{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#  Variational Inference"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The autoreload extension is already loaded. To reload it, use:\n",
      "  %reload_ext autoreload\n"
     ]
    }
   ],
   "source": [
    "import sys\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "import torch\n",
    "import torch.utils.data\n",
    "from torch import nn\n",
    "from torch.nn import functional as F\n",
    "\n",
    "import seaborn as sns\n",
    "\n",
    "from swag import data, models, utils, losses\n",
    "from swag.posteriors import SWAG\n",
    "\n",
    "import tqdm\n",
    "\n",
    "import os\n",
    "import math\n",
    "\n",
    "torch.backends.cudnn.benchmark = True\n",
    "torch.manual_seed(1)\n",
    "torch.cuda.manual_seed(1)\n",
    "np.random.seed(1)\n",
    "\n",
    "import hamiltorch\n",
    "\n",
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "metadata": {},
   "outputs": [],
   "source": [
    "def featurize(x):\n",
    "    return torch.cat([x[:, None], x[:, None]**2], dim=1)\n",
    "\n",
    "class RegNet(nn.Sequential):\n",
    "    def __init__(self, dimensions, input_dim=1, output_dim=1, apply_var=True):\n",
    "        super(RegNet, self).__init__()\n",
    "        self.dimensions = [input_dim, *dimensions, output_dim]        \n",
    "        for i in range(len(self.dimensions) - 1):\n",
    "            self.add_module('linear%d' % i, torch.nn.Linear(self.dimensions[i], self.dimensions[i + 1]))\n",
    "            if i < len(self.dimensions) - 2:\n",
    "                self.add_module('relu%d' % i, torch.nn.ReLU())\n",
    "\n",
    "        if output_dim == 2:\n",
    "            self.add_module('var_split', SplitDim(correction=apply_var))\n",
    "\n",
    "#     def forward(self, x):\n",
    "#         return super().forward(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fbd0ad7cfd0>]"
      ]
     },
     "execution_count": 86,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWoAAAD4CAYAAADFAawfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deXxU1dkH8N+TsCWAgklcAJNQS8UNRQIuVFFAxCq7UDGgRRQV2uL7qkhFxaKIVVwoi34oCAJRRFkEF6gLIr6oEKhFQVFEEkFkVREhIWSe948zkUmYfe6de2fm9/185gOZ5d4nM5lnzpzznHNEVUFERO6V5nQAREQUHBM1EZHLMVETEbkcEzURkcsxURMRuVwtOw6anZ2t+fn5dhyaiCgprV27do+q5vi7zZZEnZ+fj+LiYjsOTUSUlESkJNBt7PogInI5JmoiIpdjoiYicjkmaiIil2OiJiJyOSZqpxQVAfn5QFqa+beoyOmIiMilbCnPoxCKioAhQ4CDB83PJSXmZwAoLHQuLiJyJbaonTBq1NEkXeXgQXM9EVENTNROKC2N7HoiSmlM1E7IzY3seiJKaUzUThg7FsjMrH5dZqa5noioBiZqJxQWAlOnAnl5gIj5d+pUDiQSkV+s+nBKYSETMxGFhS1qIiKXY6ImInI5JmoiIpdjoiYicjkmaiIil2OiJiJyOSZqIiKXY6ImInI5JmoiIpdjoiYicjkmaiIil2OiJiJyOSZqShzcZ5JSFFfPo8TAfSYphbFFTe5Us/U8fDj3maSUxRY1uY+/1nMg3GeSUgBb1OQ+/nZpD4T7TFIKYKIm9wm3lcx9JilFhEzUIlJPRFaLyH9FZIOI/D0egVEKC9RKzsriPpOUksJpUZcD6Kiq5wI4D0BXEbnQ3rAopQXapX3CBGDrVsDjMf8ySVOKCJmo1Tjg/bG296K2RkWpjbu0E1UjqqFzroikA1gL4LcAJqvqPX7uMwTAEADIzc1tUxJspJ6IiKoRkbWqWuDvtrAGE1W1UlXPA9AMQDsROdvPfaaqaoGqFuTk5MQWMRER/Sqiqg9V/RHAcgBd7QmHiIhqCqfqI0dEGnn/nwHgCgBf2B0YEREZ4bSoTwGwXETWA1gD4C1Vfc3esIhchItBkcPCqfpYr6qtVbWVqp6tqmPiEZjt+OajcFRNZy8pAVSPLgbFvxeKo9ScmWjFmy+eiZ4fKs7xN52di0FRnIVVnhepgoICLS4utvy4lsnP97/QT16emUgRSs1FgwAzIcOOWt94nouOlZZmPsxrEjETb4gsEqw8LzUTdaxvvlgTfSTieS46Fp9/ipOY66iTTqC1JMJdiS3QokF2LLkZz3Mlq1i6jgJNZ+diUBRHqZmoY33zxZroIxHPcznFzj74WMcjOJ2d3EBVLb+0adNGXW/OHNW8PFUR8++cOZE9NjNT1bz1zUXE/BvpsaI5V2amtedwkt2/X15e9WNXXfLyrDk+kUUAFGuAnJq6iTpWVYneN0nblUhj+VBxOysSabDnp+Zr4/vBSuQiwRJ1ag4mWomDTbGJdWA3VFUMXx9KEBxMtBMH+2ITax98qDpnDgZSEmCijlUqDPbZKdZEGuqDkoOBlASYqH1FU33AFltsYk2k4XxQFhbauzMMZ46S3QJ1XsdyScjBxFiqD5J5sM/tnK6Kcfr8lDTAqo8gfKs3WMblrGg/8MJ9nB0fqCz/I4swUQfirzUUSxkXW9bRC9Ey3bhRddYs1RdeUH3lFdXly1X37bPo+LG8bt7yv31opN/jRJb/UdSYqAMJ1pL21zIK9obmV+DYBHgt9jQ7V4d23qRpOOL35cnPV+3fX3XxYtXDhyM/vmZlRfe6+XwT8wDaBms0A7/oUxiuR5DGFjVFjIk6kECTIfy9YUMl4mBfgdnSDs3Pa/E8Bmpj7NV0VOif8U/diJb6OU7X/9Zrp0tHvKPjxqn266eanW0ekp2tescdqt99F97xw/6ArqnG38IidFdA9Ux8poDqxWmr9IvHXrXrmaIkxUQdSKi+ad+EGqovMlgiYEs7NJ/n1wPoQxilgOqlWKGf4qygifTwYdUlS0zSrlXLPL333qv6ww/+jx/WJVjXRY1Yz8M6bYFNWoF0nZV1hzauX6YZGapTp6p6PPY8XZR8mKgDiaS7ItRU5ECJID098hZbKvK+FpUQHY6nFFAdmF6kh1ErokS6ebPpCqnq1Zg505ssA73WWVmRvz4+fwsL0UMB1VkY8GtM332n2rmzucu110bYl04pi4k6mHC7JUK1qAMlgmhabCnKM3uODqr/kgKqdzScppWzglTkhPigW7dOtX17c9fOnVW//lr9v9bRjC349E37tqZ9Y6qsVH30UdPCb97cfIAQBcNEbYVw3tD+EgHLt8L2wAPmqXngAZ8ug0gSqe/zn5WllSdk6zO4TRvKfq1f97B++WWAE0c6huCNaQF6Hm1NB4hp1SrTaD/lFNVPPw3/uaDUEyxRu2pRpilTgCNHLA8nYh4P8MMPwN69wL59ZsJcnTpA3dKv0O6/U9Fr73Q0zjvOzD4MNcuNW2mFZfZs4IYbgEGDgOnTzXP+q6Iis3ZHaamZcejveff3PHttRR5OxyYM7lSCKW//zpJ4PbOL0Hrw+ThUkY6NuVeh1iNjAr6eGzYAV1wBlJcDS5cCbdtaEgIlmWCLMrmqRR2qpDnel8aNVU87TfW3v1XNzVU94QRzfe3aqtdcY+p6w+p/ZNVHUCtWmOf08stVy8ujPEiIwcJBmK6Z8ovu3WtNzLNnm0O/+GJ49//6a1NK2KCBqQEnqgmJ0qLet8/yUKIiAjRsCNSqVf16VWDtWmDuXGDePODbb819OnYEBg4E/vhHoHZtZ2JOVN9/D5xzDpCVBXz4IdC4cZQHCrRcqtenOBut8CnGjQNGjozyHF7l5UDLlibW4mJz6nBs3w506QJ8/bX5++nePbY4KLkkTIs6kVRWqn78seo995hWN6DarJnq44+r/vST09ElBo9HtWtX1Xr1zMzDmIRRftep3kpt2jTExJgw/POf5pBLl0b+2D17VNu2NcVAs2bFFgclF3Aw0V6Vlaqvv26+ugOmq6S01Omo3G/iRPN8TZpkwcFCLQeQmalL7lyugJmGHq39+1VzcsxrHW2N9P79R/9WJkyIPhZKLkzUcbR8uepxx5n+yC1bnI7GvTZsMC3pq66ycFJIjaoPzcqqNi5QWan6u9+ZFm2053zwQfOu+fjj2EI9dEi1Z089tsqFUhYTdZytWWMGIps1U920yelo3KeiQrV1azPle8eO+J578mTzV//uu5E/9ocfzIdwr17WxFJRoXrTTSaeoUPNNzNKXcESNTcOsEFBAfDee2bQ6corgZ9+cjoid5k0CfjPf4BnnwVOPjm+5x40CDj1VOCuu8LbktHXlCnA/v3A/fdbE0utWsC0acDdd5tjX389cPiwNcem5MJEbZNWnxZhcXovfLv1CIY1XchdP7x27ABGjwa6dgV6947/+TMygHHjgHXrIntJDh4Enn7axN26tXXxiACPPWYuL70EXHMNcOCAdcenJBGoqR3LJdW7PnwHtsbgPgVU59QZxPppVR0wQLVOHQ08SzAOKitVCwpUmzZV/eWX8B5TNfC5YoV9cT33nKkGadtWLav3psQB9lF7xWviiU+pWAXStT1W6nH4Ub9p2t6e8yWIFSvM0zJqlNORqL7/vollzJjQ9z182Ex4uvhi+wf9Xn3VfJAVFKj++KO95yJ3YaJWjXzxHQt2/ai6bEG+NsRPejE+iLmGN1FVVKiec45JeOG2Yu3Wu7dq/fqq27YFv9/MmealfO21+MS1ZIlZzOmii0wpH6UGJmrVyBZHinW3Fj/negHXKaB6110W/k4JZMYM81S8/LLTkRy1ebNqRoZqhw7mg8Sf8nLV009XbdUqviV0r7xiukE6dHDPBxvZi4laNfR60r5iXfEuQKIf2nmTAqoLFlj5i7lfWZl56goKXFQv7P3G9DxuUED17qs3+L1bVd30kiXhH9OqrrUXXjCHuvZalu6lAiZq1ciSbyRJPRA/b9qyMjNQdNxxqbU+8aRJ5ulbtszpSLxqfJDehikKqM4fvqLaff57ypVaC4f1+syFYS99GvW3sADGj1fX9OuTvZioVSN7I9m4hvTWrWYyzLnnqh44EPPhXO+XX1RPPln10kstbk3H0nqt8fqWoY62xcfaUPbrm2+aDQwqMhpqG6zRHOzU3cgKezMBq/9mPB7Vm282h+LaIMmNibpKuG9um3cUf/NN1bQ01T59kv8r7T/+YZ6+lSstPGisr4+fb0wlOFVzsVUB1Va1N+oAzFJAdR6uDS/pWvEtLIDDh1U7djTVIJY+j+QqMSVqAKcCWA5gI4ANAIaHeoxrE3UkbC7le/JJ8+zff7+lh3WVn34y3x6uusriA8faeg3w+LLcFvrcc0d3E++F+eoJN+navJPPvn1mnZLsbO+2YpR0Yk3UpwA43/v/hgC+BHBmsMckRaK2mcejOniweQXmznU6Gns8/rj5/VavtvjAsbZeQ7TIPbl5ugoX6gHUuE+wpGvztzBVM0mocWPVM85gjXUysrTrA8CrAK4Idh8m6vCUl6tecolZRS7ZvtKWl6s2aWK+slvOitZrsG9M0SbdOEyoWr7c1Fh36RK4pJASk2WJGkA+gFIAx/m5bQiAYgDFubm5cf0FE9nu3aZOt1Gj5Nr8tKpuOprF9UOKQ+vVzdunTZtmfuVhw5yOhKxkSaIG0ADAWgC9Q92XLerIbN1qWp9Nmpj/J7rKStUzzzSVLbbVTbs4kcbDXXeZd+/EiU5HQlYJlqjDWj1PRGoDmA+gSFUXRLP4EwWWlwcsW2ZWaLvySmDnTqcjis3rrwMbNwIjRtTYTdxKhYXA1q1mrdKtW1NuR/dHHzV7Lg4fbv52KLmFTNQiIgCmA/hcVZ+0P6TUdPbZwJIlZsPcjh2BXbucjih6jz1mPnz69XM6kuSVnm6WaT3nHPM8b9zodERkp3Ba1O0BDATQUUQ+8V7+YHNc7lJUBOTnm+2m8/NtW1v69783rdGtWxM3Wa9eDXzwAfC//3vsLu5krQYNgMWLzRrb3boBe/Y4HRHZJWSiVtUPVFVUtZWqnue9vBGP4FyhqAgYMgQoKTHDViUl5mebkvVll5lkvWUL0KmTWWg/kUyeDDRsaHZSIfvl5gKvvgps3242YuAOMcmJO7yEMmqU6Tz2dfCgud4mVcn6m2+A9u2Br76y7VSW2rPH7FJyww0mWVN8XHABMGMGsHIlcNttpj1ByYWJOpTS0siut8jllwPLl5ttmS6+GFizxtbTWWLGDLNP5O23Ox1J6unfH3jgAfMajB/vdDRkNSbqUHJzI7veQm3bAv/3f6Z1evnlwKJFtp8yah4P8MwzQIcOwFlnOR2NTeI0VhGtBx80A4v33GP6ril5MFGHMnYskJlZ/brMTHN9HLRoAaxaZZJfr15mY9hId8+Oh2XLTFfN0KFOR2KTOI9VREPEtKjbtDE7mq9f73REZJlABdaxXJJuwosLJlccOqR6001mkkO3bmaRHje5+mqznGl5udOR2MTmRZestH272bg3N1f1+++djobChVgnvKQ8qyZXxPDVuV49YNo0YNIk4M03gVatgHffjS4Mq33zDfDGG8AttwB16jgdjU0cGquIRpMmputj926gZ0+grMzpiChWTNTxYsFXZxFg2DDTFZKZacr37rzT+TfiI48AtWsDt97qbBy2cnCsIhrnnw/Mng189JH5AGUlSGJjoo4XC8v82rYF1q0zpVhPPmlmNb7hUGX75s2mX/S224CmTZ2JIS4cHquIRp8+wEMPAXPmmCnnlLiYqOPF4q/O9eubKou33jIzAK++2nzNjXfN9Zgxprvjb3+L73njrrAQmDrVzI0XMf9Oner6NUZGjTKle/feCyxc6HQ0FC0m6nix+quzt7+7c5c0rD/UAo/+8T946y3gjDOAP/3JtHTt9vnnJoxhw4CTT7b/fI5LwIWgRIDp082kmAEDgE8+cToiigYTdbxY+dW5Rn93ndLNuGfJ7/H1Y/MxfLiZHdiyJdC3L/D22/aV8z34oPkVRoyw5/hkjYwM05o+4QSz4t733zsdEUUsUDlILJekK8+zSqAyv0jL/0KUiu3YYdYrPuEEc/VvfmP2Zly9OvbNdD0es8HBmDHm2KNGxXY8ip9168z+ChdeaMo9yV3AXchdLJrdSsLcM/DQIdWiItXLLjO7ngOm1vn661WfeEL1vfdUd+4MnLw9HrMDzYcfmgXqb7xRtXnzo6e79FL31XPHxAX18nabP9+8doWFNm7qQFEJlqhFbajbKSgo0OLiYsuPm5Ty800XRk15eaYf1KLH7NkDLF1q1rxetQrYtu3obbVqASedBDRqdLSM6+BB4Lvvqq/GduKJwIUXmoHLa64x9bpJo6o7ybcyJzMzIQYMI/Xww8D99wPjxgEjRzodDVURkbWqWuD3NiZqh6Wl+S9yFQncuWxBUtm505T4bd5sllLdsQPYv9+cVsRMsGnSxJTc5eWZaclNm9q4Y4vTovnATFCqZor5Sy+ZvusePZyOiAAmaneLNkEUFZnaq9JSUzkydmzStfziKpoPzAR26JBZQGvjRrPw17nnOh0RBUvUrPpwWrTVIAlYKuZqCTbzMFYZGWY1xkaNTCVIIu4mlEqYqJ3mO5ECMJvhVc1YdNHKbEkvAWcexqpJE7M7zO7dZneY8nKnI6JAmKjdoLDwaKKorDTXuXAZzaSWoDMPY9WmDTBzpun+uP32KNcEcfk63cmAfdRukUKDWeQ+o0eb5QDGjzcLfYUthapl7MbBxESQYoNZ5C4ej9kdZuFCU8L5hz+E+UA2MCzDwcREkGKDWeQuaWnA88+bdc779zfruIQlgdbpTmRM1G4RzmAW+wLJRvXrm8HFjAygWzdg794wHsQGRlwwUbtFqMGsBNizjxJfbq7p/vj2W9MVUlER4gEpWC3jiEBzy2O5cK0PGyTQnn2U+J5/3vx5DR2qoddASYE1UuIBXOsjCXCwkeJsxAjg8ceBZ2r/FbdVTDx6A6s6bMHBxGTAvkCKs3HjgKsz3sFfKp7Ae+hw9IZAW8hxDMU2TNSJgn2BFGfp6cALh3qjBb5CH8zHFjQ/emPNqg6OodiKiTpRpOjMOXLWcXmNsQTdoBB0x2LsR0NzQ81vchZu3kzHYqJOJFyIieJt7Ficlvk9XsG1+AItMQBzUJnR4NhvcqynthUTNREF5v0m1zFvCybgDixBd9zX+aNjGwkcQ7EVEzURBef9JjfUMwlDhgCPLjnr2K5njqHYiomaiMIiAkycCFx6KXDzzcCaNT43cgzFVqyjJqKI7N4NtGtn9tNcsybJ9s50EOuoicgyOTlmTZCffgJ69QLKyvzciTXVlmKiJqKItWoFzJ4NrF5tyqWrfTFnTbXlmKiJKCq9epnNBmbPBp580ucG1lRbLmSiFpHnRGSXiHwWj4CIyKV8uzOys4HsbNz3QBr6Zr6GEXd78Oab3vsFqp0uKWF3SJTCaVHPBNDV5jiIyM1qdmfs3Qvs3QuBYsbBP6IVPkX/aw9j0yYErp0Wqd4dMnCguY5JO6SQiVpV3wewLw6xEJFb+evO8KqPg1ik3VGnbD+6dwd+vPexY2uqRY5d/bHqZ/Zhh2RZH7WIDBGRYhEp3r17t1WHTV0cNSc3CTEVPA+lmO/pjW++Afov7IfKZ/9VvaY6VBkw+7CDsixRq+pUVS1Q1YKcnByrDpuaOGpObhPGVPBL8koxaRKwdCnwt0+vr74uTV5e6HNwXZCAWPXhRhw1J7fxN0Xcl3e6+JAhwNChZsOB2bMjeDzAdUGCYKJ2I65ERm5Tc4p4Vpa5+Jku/vTTwGWXAbfc4jPN3PfxgHmcL64LElTIKeQi8iKAywBkA9gJYLSqTg/2GE4hj1F+vunuqCkvz3yNJHK5PXuAtm3NNPPiYuCUU2rcoajIfEMsLTUt6bFjU35dkGBTyLnWhxtV9VH7dn9wnzpKMOvXAxddZGYxLl8O1KvndETuxrU+Eg1XIqMk0KoVMGsW8NFHwO3ZL0PFpgqmFKiQquV0ABRAYSETMyW8PmVFeKBWCcb8ci9aYyX+WjLRfFsErPn7rvnts6pCyqrjuwS7PojIPvn58JSUojcW4DVcg2W4Ep3wrnXjLUk0nsOuDyJyRmkp0qCYjYFoiS/QD/PMbuZWVTClSIUUEzUR2cdbG90QB/Aqevy6m/nPzc6w9PhhX5+gmKiJyD4+E11OwxbMQz98jjNw40lvwuOx9vi/SsKabCZqIrJPjQqmznmbMb7wEywszsXDD9e4bzTVGylSIcXBRCKKK1XgT38ypXsLFpgNCDh3gBNeiMhlysqADh2AjRuBDz8Ezr4mP2mqN6LFqg8icpV69UxrukEDoEcPYF/Jz/7vmGTVG9FioiYiRzRtapL1tm3AdfUW4gjSj71TklVvRIuJmogcc9FFwJQpwFtll2JkrSeq35iE1RvRYqImIkcNHgwMGwY8cWQ4irL+mtTVG9HiWh9E5LinngI++wy4+eMJaLlmAtq0cToid2GLmogcV7s28PLLwIknAj17Ajt3Oh2RuzBRE9GxHFg6NCcHWLQI2LsXuPZas+kAGUzURFSdg5srt24NTJ8OfPAB8D//E8OBkmyNaiZqIqrO4c2V+/cH7r7bVINMmxbFARz8oLELEzURVeeCpUPHjQO6dDE7mq9aFeGDA33QDBiQsK1rJmoiqs4FS4emzy3CixvPRW7FZvS5dBe+m7Qg/AcH+0BJ0NY1EzURVef00qHerosTtq3HIvTEz5WZ6DO8Kcpnvhje40N9oMSxG8cqTNREVF28lw6tOfA3fPivXRdnYwNm4QZ85LkAw/4iCGsNOX8fNDUl2BoiXD2PiJzjb3lTP+7HGDyM+zF5sum3Duu4o0b5X5EPcOWqfFw9j4jcyd/Anx9/x2hck/E2hg8HVq4M47iFhSYRz5mTFDvAMFETkXPC7IJIy8zAnAn78JvfmMkw27aFefwk2QGGiZqInBNo4C8r65jkevwt/bBoEXDoENC7t9l8ICxVrWuPx/ybYEkaYKImIicFqjCZMMFvcj3jDNObsWYNcHv2y1BJjpmHoTBRE5Fzouia6P5zEUbXGouZv/TFZAxN2NroSLDqg4gSS34+PCWl6IWFeB1X4x10Qge878pKjkiw6oOIkkdpKdKgmI2B+C02oy9eRilOTbja6EgwURNRYvEOQB6Hn7EIPVGGeuiNBTjUrIXDgdmHiZqIEovPAGRLbMIcDMBaFOC2jOeheflJs7SpLyZqIkosNQYgu+etx4MFr2HWlxdiUmm3pFna1BcHE4ko4XnymqNX6dPVBxeBhBpg5GAiESW1tG9LMAs3/Dq4+C2amRuSZICRiZqIEl9uLo7H/l8HF/tgPspQN65raNuJiZqIEp93gLElNmE2BmIN2mFo+lTow4m1+FIgYSVqEekqIptEZLOIjLQ7KCKiiPgMMPaQJbj/+AmYUXkDnv058db18CfkYKKIpAP4EsAVALYBWAOgv6puDPQYDiYSkZM8HqBbN+Df/wbeew9o397piEKLdTCxHYDNqrpFVQ8DmAugh5UBEhFZKS3NVOY1b26WRf3uO6cjik04ibopgG99ft7mva4aERkiIsUiUrx7926r4iMiikqjRsDChcDPP5tkffiw0xFFz7LBRFWdqqoFqlqQk5Nj1WGJiKJ21lnAjBnAhx+arRgTVTiJejuAU31+bua9jojI9fr2Be65B3j2WWD6dKejiU44iXoNgBYi0lxE6gC4DsBie8MiIrLO2LHAFVeYjXFXr3Y6msiFTNSqegTAnwEsA/A5gHmqusHuwIiIrJKeDrz4ItCkCdCnD7Brl9MRRSasPmpVfUNVf6eqp6lqclSQE1FKycoCFiwA9uwB+vUDKiqcjih8nJlIRO5WVGSWLbVg+dLWrYF//QtYsQIYMcKyCG1Xy+kAiIgCKioyy5UePGh+rlq+FIh6N/EBA8zmuE8/DbRtC1x/vUWx2ojLnBKRe+Xnm+RcU4zLl1ZUAJ06AcXFwKpVwHnnRX0oy3CZUyJKTIGWKY1x+dLatYF584DGjYHevYF9+2I6nO2YqInIvQItU2rB8qUnnwzMnw9s3266PyorYz6kbZioici9fPZH/FVmprneAhdeCEycCCxbBowebckhbcFETUTuVWN/ROTlmZ+jHEj055ZbgMGDTe5ftMiyw1qKg4lElPLKyoBLLgE2bTIVIaefHv8YOJhIRBREvXqmv7puXaBXL7PinpswURMRwYxPvvSSaVXfdBNgQ2dD1JioiSgxWDhDMZCOHYF//AN45RVg/HjLDx81zkwkIvezYYZiIHfeaVbYGzkSOP98MzHGaRxMJCL3s2mGYiAHDgAXXGBW2Vu71pKy7ZA4mEhEic2mGYqBNGhgVtorLzfLopaV2XKasDFRE5H72ThDMZDTTwdmzTLrgfzlL7adJixM1ETkfjbPUAykZ0/g3nuBadPMxSlM1ETkfnGYoRjImDFAly7AsGFmMowTOJhIRBTC3r1AmzaAxwOsWwdkZ1t/Dg4mEhHFoGobr127gP7947/SHhM1EVEYzj8feOYZ4O23gfvui++5maiJiMI0aBBw663Ao48CCxfG77xM1EREEZgwAWjXDrjxRrMuSDwwURMRRaBuXbMWSN26ZhuvAwfsPycTNRFRhE49FZg7F/jiC7PpgN0r7TFRExFFoVMn4JFHzCa5EybYey4maiKiKI0YYTYauOsuYOVK+87DRE1EFCURYMYM4LTTgH79gB077DkPEzURUQyOP95Mhtm/H+jbF6iosP4cTNRERDE66yxg+nTgzDPNNHOrcYcXIiILXHedudiBLWoiIpdjoiYicjkmaiIil2OiJiJyOSZqIiKXY6ImInI5JmoiIpdjoiYicjlbNrcVkd0ASiw/sHWyAexxOggHpfLvz989dbn9989T1Rx/N9iSqN1ORIoD7fabClL59+fvnpq/O5DYvz+7PoiIXI6JmojI5VI1UU91OgCHpfLvz989dSXs75+SfdRERIkkVVvUREQJg4maiMjlUipRi0hfEdkgIh4RKahx299EZLOIbJctJ+0AAAKFSURBVBKRK52KMR5E5EER2S4in3gvf3A6JruJSFfva7tZREY6HU+8ichWEfnU+3oXOx2P3UTkORHZJSKf+Vx3goi8JSJfef9t7GSMkUipRA3gMwC9Abzve6WInAngOgBnAegKYIqIpMc/vLh6SlXP817ecDoYO3lfy8kArgJwJoD+3tc81Vzufb0TspY4QjNh3su+RgJ4R1VbAHjH+3NCSKlEraqfq+omPzf1ADBXVctV9RsAmwG0i290ZKN2ADar6hZVPQxgLsxrTklKVd8HsK/G1T0APO/9//MAesY1qBikVKIOoimAb31+3ua9Lpn9WUTWe78iJsxXwCil4utbkwL4t4isFZEhTgfjkJNUdYf3/98DOMnJYCKRdJvbisjbAE72c9MoVX013vE4JdjzAOAZAA/BvHkfAvAEgJviFx054Pequl1ETgTwloh84W11piRVVRFJmNrkpEvUqto5iodtB3Cqz8/NvNclrHCfBxH5F4DXbA7HaUn3+kZKVbd7/90lIgthuoNSLVHvFJFTVHWHiJwCYJfTAYWLXR/GYgDXiUhdEWkOoAWA1Q7HZBvvH2mVXjCDrMlsDYAWItJcROrADBwvdjimuBGR+iLSsOr/ALog+V9zfxYDuNH7/xsBJMw37KRrUQcjIr0ATASQA+B1EflEVa9U1Q0iMg/ARgBHAAxT1UonY7XZYyJyHkzXx1YAtzobjr1U9YiI/BnAMgDpAJ5T1Q0OhxVPJwFYKCKAec+/oKpLnQ3JXiLyIoDLAGSLyDYAowE8CmCeiAyGWYa5n3MRRoZTyImIXI5dH0RELsdETUTkckzUREQux0RNRORyTNRERC7HRE1E5HJM1ERELvf/ZSNx7OSJVeYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "arr = np.load(\"data.npz\") \n",
    "x = torch.from_numpy(arr['x'])\n",
    "f = featurize(x)\n",
    "y = torch.from_numpy(arr['y']) * 10\n",
    "\n",
    "x_ = torch.from_numpy(arr['x_'])\n",
    "f_ = featurize(x_)\n",
    "y_ = torch.from_numpy(arr['y_']) * 10\n",
    "\n",
    "plt.plot(x.data.numpy(), y.data.numpy(), \"ro\")\n",
    "plt.plot(x_.data.numpy(), y_.data.numpy(), \"-b\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_steps = 3000\n",
    "\n",
    "net = RegNet(dimensions=[10, 10, 10], input_dim=2)\n",
    "optimizer = torch.optim.SGD(net.parameters(), lr=5e-6)\n",
    "criterion = torch.nn.functional.mse_loss\n",
    "noise_var = 0.0005\n",
    "prior_var = 100\n",
    "\n",
    "for epoch in range(n_steps):\n",
    "    optimizer.zero_grad()\n",
    "    preds = net(f)\n",
    "    loss = criterion(preds, y) / (2 * noise_var)\n",
    "    loss.backward()\n",
    "    optimizer.step()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 109,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fbceeed7c90>]"
      ]
     },
     "execution_count": 109,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD4CAYAAADvsV2wAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deZhU1bX38e9iCNCiNwqtELAbozhjVDqIaF4FgiJRiCCItkMSFYc4JA43JkQcueq9xinGKFGvRlpFiUacJQ4xxqA2BhAc0QiIRkC8KKJA0+v9Y1dL01b1UHWqTg2/z/Ocp6tOnT5nV1X3ql3r7LO2uTsiIlL82sXdABERyQ0FfBGREqGALyJSIhTwRURKhAK+iEiJ6BB3A1Lp3r279+nTJ+5miIgUlNmzZ69w9/Jkj+VtwO/Tpw+1tbVxN0NEpKCY2aJUjymlIyJSIhTwRURKhAK+iEiJUMAXESkRCvgiIiVCAb/Q1dRAnz7Qrl34WVMTd4tEJE/l7bBMaYWaGpgwAdasCfcXLQr3Aaqr42uXiOQl9fAL2cSJG4N9gzVrwnoRkSYyDvhm1tnMXjKzuWa2wMwuTrLNj8xsuZnNSSwnZnpcARYvbtt6ESlpUaR01gJD3H21mXUEnjezx9x9VpPtprn76REcTxpUVIQ0TrL1IiJNZNzD92B14m7HxKJptHJh8mQoK9t0XVlZWC8i0kQkOXwza29mc4BlwEx3fzHJZmPMbJ6ZTTezbaM4bsmrroYpU6CyEszCzylTdMJWRJKyKOe0NbNvAg8AZ7j7/EbruwGr3X2tmZ0MHOnuQ5L8/gRgAkBFRUX/RcnSFSIikpKZzXb3qmSPRTpKx93/D3gGGN5k/cfuvjZx9xagf4rfn+LuVe5eVV6etLqniIikKYpROuWJnj1m1gUYBrzRZJueje6OBF7P9LgiItI2UYzS6QncYWbtCR8g97r7w2Z2CVDr7jOAM81sJFAHrAR+FMFxRUSkDSLN4UepqqrKNQGKiEjb5CyHLyIi+UsBX0SkRCjgi4iUCAV8EZESoYAvIlIiFPBFREqEAr6ISIlQwBcRKREK+CIiJUIBX0SkRCjgS+mpqYE+faBdu/CzpibuFonkRBTF00QKR00NTJiwcfL3RYvCfdDEMVL01MOX4ta0N3/WWRuDfYM1a2DixDhaJ5JTRRfw3eGaa+Cjj+JuicSuoTe/aFH4w1i0CD7+OPm2ixfntm0iMSi6gP/223D++bDTTnDTTVBfH3eLJDYTJ369N59KRUV22yKSB4ou4O+4I8ybB/37w6mnwqBBMGdO3K2SWLS2115WBpMnZ7ctInkgiikOO5vZS2Y218wWmNnFSbbpZGbTzGyhmb1oZn0yPW5zdtoJ/vIXuPNOePfdEPx//nP47LNsHlXyTqpee7duUFkJZuHnlCk6YSslIYoe/lpgiLt/B9gTGG5mA5tscwLwibvvAFwDXBnBcZtlBsccA2++CSedBNddBzvvDNOnh3SulIDJk0PvvbGysvDH8N57Id/33nsK9lIyMg74HqxO3O2YWJqG1FHAHYnb04GhZmaZHrs1ttwy5PJfeAHKy2HsWBgxAt55JxdHl1hVV4feu3rzIkBEOXwza29mc4BlwEx3f7HJJr2AJQDuXgesArol2c8EM6s1s9rly5dH0bSvDBwItbVw7bXw97/D7rvDpZfC2rWRHkbyTXW1evMiCZEEfHff4O57Ar2BAWa2e5r7meLuVe5eVV5eHkXTNtGhQxiG/frrMHIkTJoE/frBzJmRH0pEJO9EOkrH3f8PeAYY3uShpcC2AGbWAfgPIMWA6Ozr1QumTYMnngj5/IMOgiOPhKVL42qRiEj2RTFKp9zMvpm43QUYBrzRZLMZwPGJ20cAT7vHf+r0oIPg1Vfh4othxoxwUvfqq2H9+rhbJiISvSh6+D2BZ8xsHvAyIYf/sJldYmYjE9vcCnQzs4XA2cD5ERw3Ep07h9TOggXw//4fnHMODBkCX3wRd8tERKIVxSidee6+l7vv4e67u/slifWT3H1G4vaX7j7W3Xdw9wHu/m6mx43at78NDz8Md9wRTupWV8OGDXG3SoqKqnRKzIruSttMmMFxx4VaPA88EC7WajHxVCj/xIXSzmKVrK7PhAl6HyS33D0vl/79+3uczj7bHdyvuqqZjaZOdS8rCxs2LGVlYX1rTZ3qXlnpbhZ+tuV323KMTNspmams3PT1b1gqK+NumRQZoNZTxFXz+M+dJlVVVeW1tbWxHb++HsaPD1fmTp8Oo0cn2ahPn9BTa6qyMoz5bknT2uwQrgSN+uKgTNspmWvXLvnXRTNV+JNImdlsd69K+pgCfmpffBFO4M6dC88+CwMGNNkg03/iXAViBZv46UNXcqS5gK8cfjO6dIEHH4QePeCww5L8X6YqztXaUrupqjlGXZs903ZK5lLV9VGVTskhBfwWbL01PPIIrFsXavCsWNHowUz/iXMViEsh2OTipHQmx1BdH8kHqZL7cS9xn7Rt6tln3Tt3dt9rL/eVKxs9kMlJ11yeTM3FyeG45OJ11IlvKRDopG00Hn8cRo2CvfYK9Xc23zyCndbUhJmZFi2C9u3D4P/KytD7Vu+vdXKRH1cOXgqEcvgRGT4c7r03VN38wQ9g9eqWf6dF1dUbUy4NV3ppjHbbRHUupLmUTa7Ot4hkkQJ+G40aFeLACy+EWjyrVkWw02Rzr65ZE9ZLy6I4F9LShVE68S1FQAE/DUceGapt1tbC978PK1dmuEP1HjMTxUnplj50S+HEtxQ9Bfw0jRkD998fJkwfPBg++iiDnan3mJkoRsC09KGbi1E2Kn8h2ZbqbG7cS76N0knlySfdu3Rx33FH98WLEyvbOiJGI0DiF3fpA/0NSERoZpRO7IE91VIoAd/d/W9/c9+iy1qvaL/E32b7EOjb+o9bzMMmcy2d1zLugBv3B44UDQX8bJs61Wd3HuTdWeY9+MDn0k//uHHJJHC39oMiGx/OTTsJDYtZ5vuWktJcwNc4/Cgkxmi/zs4MYyar6crDHMr+/H3jNm2pW9MwNn/x4pDH15j8Vlm3Dt7f7nss/qA9S+lFO+r5Jv/HlnzCdr3r2GZJBH9PzRW8g/TfN43zl4g0Nw4/4544Ya7aZ4DXgAXAWUm2ORBYBcxJLJNa2m9B9fAb9c7eo8J35A3vwuf+MCOS9/Cb6yHGnVooQEuWhHLWXbsm7yQ3LD3af+TDedTP3eImn3LCP/zZZ90/+qiNB0uVeunWLb33reFvoaE3r/ddMkQ2UzqEKQ73TtzeHHgL2LXJNgcCD7dlvwUV8JsEgWV09/687O1Z77dz3Kb/uC0FdOVyW23xYvef/MS9Y0f39u3djz7a/bZu5/pMhvrr7OSvsbO/wEB/mBF+DT/z4/lf34M53okvNnlpy8vdBw92P+ss99tvd587133duhQHTZV6SbU0974l+1to2L/O40iashrwv7ZDeBAY1mRdcQf8JP+4q9jChzLTwf2KI1/x+vrEti0F9OYCik7qurv7Z5+5X3BBGB3VqZP76ae7/+tfiQebC6KJpY52/i8q/fGtj/Vrr3U/8UT3ffbZ9Nc6dHDfeWf3UaPczz/f/a673BcscF9f8e22BfzmcvD6cJcsyFnAB/oAi4Etmqw/EPgYmAs8BuyW4vcnALVAbUVFRdZfmEglSdOsXet+1FHhVT7jDPe6Om/55FyqIKCv+758ufvkye49eoSXYPx49/feS7Jh0/eilcG4rs79tdfca2pCkD/8cPdddgnB/6uXvdN6P6Td434NZ/l8dvX6hveiW7e2B2+dqJUsyEnAB7oCs4HRSR7bAuiauD0CeLul/RVUD78ZGzZsnC5x9Gj3Ndvu2HxgaEUPtdR6gq+8ElI3nTqFpz1smPsLL7RhBxn2pNeuDWmeP/7R/ac/dd+xx6qvdtGt3cc+cu/FfnV1ra/u0r1tH8rq4UsWZD3gAx2BJ4CzW7n9e0D35rYploDf4JprQtzet+8yX9alovnAkGYPtZisWxdehn33DU+1Sxf3k08OaZU2y8JQzUWL3G+7LXwQ9e0bdrld+af+9DbjdcGdxCrbJ20N+CNwbTPb9GDjdIoDEmkfa26/xRbw3d2nTw819XfYZpW/9a0DWh8YSqgn+OWX7r///canvMMO7ldf3WQOgnS0Zex8GiNnnnsutBXcTzmlDe3VBXcSsWwH/P0BB+Y1GnY5AjgFOCWxzemJIZtzgVnAoJb2W4wB3z2kIrp3d99qqzCpSquUQE+wrs791lvde/UKT2+ffdwfeiikxHIq2Wvdyg/azz8P6Tuz8P5ee21IB4nkUk5H6US1FGvAd3dfuDCMAOnYMQwDbJUi7gn+9a9hJjEIKZyZM33jqKZcay6F1spU2j//6T50aNi0b1/3G25wX7EiN80XUcDPQytXug8ZEt6BX/4yhp5sHli2zL26OrwG227rfvfdMQb6Bq0ZZ9+KVFp9vfsjj7jvvXf4lW98w33MGPd587L/FKS0NRfwVR45JltuGaZMnDABLr8cRo+Gzz6Lu1W54Q5Tp8Iuu4QZxC64AN54A8aPDxUoYtVSSepW1sA3C5Pez54Nc+bAaafBM89A//7h/a6ri6i9Im2R6pMg7qXYe/gN6uvdr7/evV079913d3/33bhblF3Llrkfdljo9Q4c6D5/ftwtaiKLV78uX+4+blzY1YABIbUnEjXUw89fZnDGGaG3//77UFUFTz0Vd6uy4/HHoV8/ePJJuPpqeP552G23uFvVRLKJTu68M4T+997LqIhd9+5hprRp02DhQvje9+D116NrukhLFPDzxLBh8NJL0KNHmCv3mmtCjCkG69fDOefAIYdAeTm8/DL8/OfQvn3cLUuhujoE9/r6jIN8MuPGwXPPhd0feCDMnx/p7kVSUsDPI337wqxZYaL0s8+GY4+Fzz+Pu1WZ+fBDGDo09Oh/+tMQ7Pv1i7tV8dttN/jrX6FDhxD058yJu0VSChTw88zmm8P06XDppXDXXbDvLp+wsNcBBTnP6d//DnvvHU5c1tTADTdA585xtyp/7LRT6OmXlYVvdW+9FXeLpNgp4Oehdu3g17+Gx857mqVL6qn64EEe9MPCBBkTJhRE0L/zzjC5e9eu4VvL0UfH3aL8tP328Je/hNvDhoXzOCLZooCfxw6e9hNm058dWMgPeZDz+G/Wr1kXZlXKU/X1YZjlccfB/vvDiy8qhdOSHXcMJ7Q/+ST09FesiLtFUqwU8NNRUxPSK61Js7Rl26YWL6YPi3ie/TmN33EV53Egz7JkUSunSsyxL78M5zcvuwxOOCEEsa22irtVhWHvveGhh+Ddd2HkyPBaikQu1XjNuJe8HYfflro2mdbAaXKZ/90c6V351Ldqt9JnzIj0WWVsxQr3/fYLTb3iijy4YrZA3XdfeA2rq/UaSnpQaYUItaVyZaZVLpN8YLzZeQ/fs/JjhzAl35dfRvfU0rVwYagZ06mT+7Rpcbem8F12WXi7L7007pZIIWou4Cul01aLF7d+fVu2TSbJRUA73vKfzHpzK848E667DgYOjPfinXvvDReLrVwZLhgbNy6+tsQik5RdCr+qrOGYze7nggvgvvLTCuIkvRSIVJ8EcS/q4bfswQdDqeXOnd1vvDG3KYBVq9yPP96/KmVckmUCWpOya2uV08Q+v+QbPojnvYzVPr9z/6KqjirZhVI6EcplDr8VPvjA/aCDwq4PPdT9ww8j23VKs2e7f/vbof7PpElhdqqClUnZ6ZY+0NN5/xvtcyk9fWv+7Tvxun+67a5pP0UpLQr4UUtn9qQs1rHfsCFMttGpU5hLe/r0yA/h7uEbxM03h+P07u3+/PPZOU7OZPqBnO6E9G2Y2PwZDvB21Pk47tFJXGmVrAZ8YFvgGeA1wqxWZyXZxoDrgYWEmbH2bmm/eR3w89SCBe79+4d39ZhjIpgWsJE1a9yPOy7s+6CDQtXLgpdpyq2l32/pA6GV+7yC/3QIH+oiLcl2wO/ZEMCBzYG3gF2bbDMCeCwR+AcCL7a0XwX89Kxb537hhe4dOrj37Bny/Jl6/333qqrw13LhhWE6wqKQTkBurKVvCOl8oCTZZ32XMh+592Lv2NH9lVcyfdJS7HKa0gEeBIY1WXczcFSj+28CPZvbjwJ+Zl55xf073wnv8NFHh1rs6Zg1K3xwdO3q/uc/R9vG2EVxUr25lF26KaMk+1yxwv1b3wpTY37+eVufqJSSnAV8oA+wGNiiyfqHgf0b3X8KqEry+xOAWqC2oqIi269L0Vu71v2ii8Lcud27hzjSljzwvfeGfP1227m/+mr22hmbXEwOH+E5nJkzQxNPOSWy1kkRyknAB7oCs4HRSR5rVcBvvKiHH51XXw1DJ8F9+PCWZ9Wqr3e/6qqw/X77pf/toCAU2OTw554b3pei+7YlkWku4Edy4ZWZdQT+BNS4+/1JNllKOLnboHdineTA7ruHUsXXXbdxlqnLL4d1676+bV0dnHkmnHsujB0bKjl27577NudMlic7idpll8Gee8KJJ8Ly5XG3RgpNxgHfzAy4FXjd3a9OsdkM4DgLBgKr3P3DTI8trde+fQjkr78eJtf+1a9C4HjkkZDLAFi6FIYMCXXrzz4b7rlH9evzTadOofT0qlVw+ulxt0YKTRQ9/P2AY4EhZjYnsYwws1PM7JTENo8C7xKGZf4BOC2C40oaevcOE6w88kiYevDQQ0Pd+ilTYK+94JVXYOpU+M1vQrUAyT+77w4XXhjKWtyf7Pu0SArmDd27PFNVVeW1tbVxNyM6NTWhjv3ixVBRAZMnx54+WL8e/vAHuPhiWLYspHqmT4edd461WdIK69fDPvuEb2WvvQbdusXdIskXZjbb3auSPaY+XC7U1ISZqhYtCvmTPJm5qmNHOO00WLgQpk0Lk6gr2BeGjh3h9ttD0bozz4y7NVIoFPBzYeJEWLNm03Vr1uTNzFWbbx6qXJaVxd0SaYs99ghTYd51FzzxRNytkUKglE4utGu38cxoY2ZhdIhImtauDYG/rg7mz4cuXeJukcRNKZ24VVS0bb0Utwhr6HfqBDfdFKZG/K//iqyFUqQU8HNh8uSv50vKysJ6KS1ZOJ8zeDAceyxceWW8k+FI/lPAz4UkM1cxZUrso3QkBlk6n3PVVbDZZnDyyWEEj0gyCvi5EuUVnVmYVk9yJNNpL1PYemu49lr4299g9Gj44ouMdidFSgG/0OTpEE9ppSyezzn+eLjxxnBR3Q9+AJ99lvEupcgo4BeaPB/iKS3I8vmcU08NpReeew6GDoU33ohkt1IkFPALTZZSApIjOTifU10dSi689Rb06xcK4X36aWS7lwKmgF9ook4J6HxA7uWgQufIkfD22/DjH8PVV8OOO4YLtPL0shvJEQX8QhNlSkDnA4paeXn48vDyy+GLRHU1HHRQ+CCQ0qSAX2iiTAnofEBJ6N8fXngBfve7UC+pX78wH4KGb5YeBfxClCwlkE5qRucDsiMP02Tt24dCeW+8AYcdFuZD2GcfmDMn7pZJLingF4N0UzMq+RC9PE+T9ewJ990XymB/8AF897uhPLZ6+6VBAb8YpJuaUcmH6BVImmzMmFBH/8gj4aKLYNAgDeEsBVHNaXubmS0zs/kpHj/QzFY1mhFrUhTHlYR0UzMq+RC9AkqTbbVVmN3svvvgX/8KM5795jewYUPcLZNsiaqHfzswvIVt/ubueyaWSyI6rkBmqZkCm8Q77xVgmuyII0Jp5WHDwpj9QYNgwYK4WyXZEEnAd/fngJVR7EvSoNRM/ijQ96JHD3jwQbj77lBqea+94NJLldsvNrnM4e9rZnPN7DEz2y3ZBmY2wcxqzax2+fLlOWxagWuamunWLcyEceyxeTNKpGQUcJrMDMaPD7n9MWNg0iQYMEAjeYqKu0eyAH2A+Ske2wLomrg9Ani7pf3179/fJQ1Tp7qXlbmHMSJhKSsL60Xa4IEH3LfZxr1DB/eLLnJfty7LB5w61b2y0t0s/NTfbFqAWk8RV3PSw3f3T919deL2o0BHM+uei2OXnAIZJSL574c/3HQkz377ZXEkT54PZy0WOQn4ZtbDzCxxe0DiuB/n4tglp4BGiUj+azyS5513Qm7/hhuyUJNHHZWciGpY5t3AP4CdzOx9MzvBzE4xs1MSmxwBzDezucD1wPjEVw+JWgGOEpH81zCSZ/BgOOOMUJwt0tNs6qjkRFSjdI5y957u3tHde7v7re5+k7vflHj8Bnffzd2/4+4D3f2FKI4rSRToKBHJfz17hslVrrsOnnwS9tgDZs6MaOfqqOSErrQtNq0ZJZKHtV6kMJjBmWeGCpxbbgkHHwwXXAB1dRnuWB2V3Eh1NjfuRaN0skSjeCQiq1e7//jH4U/ogAPcly7NcIcapRMJ4h6lI3lEJ8ckIpttBrfdBnfcEXr8e+8dyjCn1NI3S131nXUK+KVGJ8ckYscdF+rsb745HHgg3HJLko3aMuxSKcesUcAvNTo5Jlmw224h6A8eDCedFPL89fWNNmjtN0uNx88qBfxSo5NjkiVbbhlG8fz85/Db34aMzLp1iQdb+81SKcesUsAvNQVc60XyX4cOYdL0K66Ae+4JV+uuWUPrv1kq5ZhVCvilSCfHJMt+8YvQj3j88TB089OJV7bum6VSjlmlgC8iWXHSSaGXP2sWDLvlSD655vaWv1kq5ZhVHeJugIgUr3HjoHNnGDsWhvx+LDNrx9K9ubKJDR8AEyeGNE5FRQj2+hYaCfM8LWlTVVXltbW1cTdDRCLwxBMhn7/99vD007D11nG3qHiZ2Wx3r0r2mFI6IpJ1Bx8Mjz4aZtMaMiTiwmvSagr4IpITgwfDww8r6MdJAV9EcmbIEHjoIVi4EIYOhY9bmhVDV91GSgFfRHJq6NAQ9N96K6R6Vq1KsaGuuo2cAr6I5FZNDd8/sQ9/Wnso82avY8R3l7N6dZLtdNVt5KKa8eo2M1tmZvNTPG5mdr2ZLTSzeWa2dxTHFZEC06jX/gMe4S6OZtbbWzFq86f4snKnTXvvqa6uXbRIaZ40RdXDvx0Y3szjhwB9E8sE4PcRHVdECkmTXvsR/Inb+RFPM5Txi6+k7qRTNwbwVFfXmm2a5jn22LBOwb9FUU1x+BywsplNRgF/TNTnnwV808x6RnFsESkgSXrtxzKV6zmDB/khJ31xHf6rRMom2VW3Zl+fQb3hvnL8LcpVDr8XsKTR/fcT6zZhZhPMrNbMapdrzFY0NMpB8kmKXvsZ3MCFXMTt/JjzFp8RYniyQn8tXSiqHH+z8uqkrbtPcfcqd68qLy+PuzmFT6McJN8k67UnXMjFnM5v+Q3ncPXViZVNC/1VVrZ8DFXWTClXAX8psG2j+70T6ySbNMpB8k3jXjuEnnuCAdd1+SVHDFjEuefC9OlJfr+ZD4yvqLJmSrkK+DOA4xKjdQYCq9z9wxwdu3Sptrjko4ZeuzvceecmKZt2f7iZPz5byaBBcMwxSebIbeYDA1BlzRZEUjzNzO4GDgS6Ax8BFwIdAdz9JjMz4AbCSJ41wI/dvdnKaCqeFoE+fUIap6nKyvAPJ5KnVqyAQYNg5cpQXnmHHVJsWFOjyppNNFc8TdUyi1lDDr9xWqesTDNcSUF45x3YZ59QWXPWLNhii7hbVBhULbNUaTpDKWDbbw/33RdKMFRXw4YNcbeo8CngFztNZygFbPBguP76UGXzgsPnZ3eIcQkMYdaMVyKS1049FeZOf5vLH9qdPRjIeBZtHGIM0XRimqY/o95/nlAOX0Ty3rrKvgxZ/L/8k72YxUD6kSjbFdUAhCIa4KAcvogUtG8seYf7GMt/sIrDeYBP+GZ4IKohxiUyhFkBX0TyX0UFPfk30zmCxVRwDFOpx6K7yCrVforsIi4FfBHJf4krbAfxD67jLB7lB1za8ZLoLrJKdgVvEV7EpYAvIvmv0RDjU7iZ4zb7ExfXTeTJ8iQnVNMZbVMiQ5h10lZECs7nn8PAgfDvf8M//wm9eyce0MWGOmkrIsVls81CcbUvv4Rx42D9+sQDKhjYLAV8ESlIO+0Et94K//gH/OxniZUlMtomXbrwSkQK1rhx8PLLcNVVsMsucHpFRfLx9EU22iZd6uGLSEG74goYORLOOgseH3dbSYy2SZcCvogUtPbtw7nafv3gyJuHMP+CaUU/2iZdCvgiUvC6doWHHgqd+e9deShP3PyeCgYmoYAvIkVh223DCdyKChgxIuT183TUeWwiCfhmNtzM3jSzhWZ2fpLHf2Rmy81sTmI5MYrjiog01qdPmBZxzBg47zw46ij49NO4W5U/Mg74ZtYe+B1wCLArcJSZ7Zpk02nuvmdiuSXT44pIAYihxvxmm8G0aXD55WGsfv/+8MorWT9sQYiihz8AWOju77r7OuAeYFQE+xWRQtZw1euiRSG30lBjPgdB3wzOPx+efRa++AL23TcM1Fm1qo07KrJJUaII+L2AJY3uv59Y19QYM5tnZtPNbNtkOzKzCWZWa2a1y5cvj6BpIhKbPLjqdf/9Yc6ckNP/9a9Dfv8Xv2jldVgxfmBlS65O2j4E9HH3PYCZwB3JNnL3Ke5e5e5V5eXlOWqaiGRFnlz12r07PPAA1NbC8OHhZG5lJey2G5x9Njz2WKjN8zWpPrCOOaZge/tRBPylQOMee+/Euq+4+8fuvjZx9xagfwTHFZF8lg815hulZPqP6cO0kTW8/Tb8z/9Ar15w442h97/VVmH+3MmTYdYsqKuj+Q+mQu3tu3tGC6E8w7vAdsA3gLnAbk226dno9uHArJb2279/fxeRAjZ1qntZmXtIiISlrCysz5Pjf/65+xNPuJ93nvuee27cbIst3A/r8qRfzc/8n3zHN2Cb7qdhqazMzXNpA6DWU8XrVA+0ZQFGAG8B7wATE+suAUYmbl8OLEh8GDwD7NzSPhXwRYrA1KkhKJqFn9kM9k2P1a1bm4P0smXu06a5T5jgvsM2q776lW4s9/Hc5XdwrH/INhv3ZZa955Om5gK+6uGLSOFLVgc/FbNwFW4rLLn+AZ695DlmfrwXT3Awy9gGo579eZ5x3MuYXi/S8/2XM2x8tJqrh6+ALyKFr0+f5FUyk6msDCUX2qKmhvqTTmbuF1/Qu8gAAAejSURBVH15iMO4j7HMpx9mzve+Z4wbFy726tGjrQ2PniZAEZHi1tqRP+lWzqyupt0fbmavyk+YZJfxauVhLLjyISZNMlasgNNPh299Cw44AG64AT78sO2HyAX18EWk8KXq4XfrFiqrLV4cRgdNnpyVYmoLFsB994XltddC1mi//WDsWBg9utEUjDmglI6IFLc8msv2tddCSYfp0+HVV8O6ffeFI44IaZ/KyuweXykdESlu1dUhuLe1Dn4WSifsuitMmgTz5sEbb8Bll4W5d885JxyiqirU+XnrrYwP1Wbq4YtIacrxt4J33oE//SksL70U1u2+e+j1jxkTbptlfhyldEREmkqV909nFE8bLVkC998fgv/zz4dB/TvsEPL9o0fDd78bvnSkQwFfRKSpdu2Sz5DShnH6UfjoI/jzn0O9n6eeCmUdBgyAF19Mb3/NBfwOmTRURKRgVVQk7+HnstYPsM02cPLJYfnkE3jkkZDzzwadtBWR0jR5csjZN1ZWFqqpxVQDf8stQzHOE7M0J6B6+CJSmhpOzE6cuHGc/ogRcMcdG0/kNlTFbLx9AVMOX0SkQYwncqOicfgiIq2RJ5O2ZIsCvohIg3yYtCWLFPBFRBqkOpGbTsG1PKSALyLSIN0SDQUikoBvZsPN7E0zW2hm5yd5vJOZTUs8/qKZ9YniuCIikauuDido6+vDzyIJ9hBBwDez9sDvgEOAXYGjzGzXJpudAHzi7jsA1wBXZnpcERFpmyh6+AOAhe7+rruvA+4BRjXZZhRwR+L2dGCoWRRlgkREpLWiCPi9gCWN7r+fWJd0G3evA1YB3ZruyMwmmFmtmdUuX748gqaJiEiDvDpp6+5T3L3K3avKy8vjbo6ISFGJIuAvBbZtdL93Yl3SbcysA/AfwMcRHFtERFopioD/MtDXzLYzs28A44EZTbaZARyfuH0E8LTna00HEZEilXHxNHevM7PTgSeA9sBt7r7AzC4Bat19BnArcKeZLQRWEj4UREQkhyKplunujwKPNlk3qdHtL4GxURxLRETSk1cnbUVEJHsU8EWkNNTUxDaxSb7QBCgiUvxqasJEJkU6sUlrqYcvIsVv4sSNwb7BmjVhfQlRwBeR4lfkE5u0lgK+iBS/Ip/YpLUU8EWk+BX5xCatpYAvIsWvyCc2aS2N0hGR0lBdXXIBvin18EVESoQCvohIiVDAF5HSUsJX3CqHLyKlo8SvuFUPX0RKR4lfcauALyKlo8SvuM0o4JvZVmY208zeTvzcMsV2G8xsTmJpOhuWiEhulPgVt5n28M8HnnL3vsBTifvJfOHueyaWkRkeU0QkPSV+xW2mAX8UcEfi9h3ADzPcn4hI9pT4FbeZjtLZxt0/TNz+N7BNiu06m1ktUAdc4e5/zvC4IiLpKeErblsM+Gb2F6BHkoc2Oa3t7m5mnmI3le6+1My+DTxtZq+6+ztJjjUBmABQUSI5NRGRXGkx4Lv791M9ZmYfmVlPd//QzHoCy1LsY2ni57tm9iywF/C1gO/uU4ApAFVVVak+PEREJA2Z5vBnAMcnbh8PPNh0AzPb0sw6JW53B/YDXsvwuCIi0kaZBvwrgGFm9jbw/cR9zKzKzG5JbLMLUGtmc4FnCDl8BXwRkRzL6KStu38MDE2yvhY4MXH7BaBfJscREZHMmXt+psrNbDmwKO52tKA7sCLuRsSklJ87lPbzL+XnDvn//CvdvTzZA3kb8AuBmdW6e1Xc7YhDKT93KO3nX8rPHQr7+auWjohIiVDAFxEpEQr4mZkSdwNiVMrPHUr7+Zfyc4cCfv7K4YuIlAj18EVESoQCvohIiVDAT4OZjTWzBWZWb2ZVTR77pZktNLM3zezguNqYC2Z2kZktbTS5zYi425RtZjY88d4uNLNU8z8ULTN7z8xeTbzftXG3J5vM7DYzW2Zm8xuta9WkT/lKAT8984HRwHONV5rZrsB4YDdgOHCjmbXPffNy6ppGk9s8GndjsinxXv4OOATYFTgq8Z6XmsGJ97sgx6K3we2E/+PGWjvpU15SwE+Du7/u7m8meWgUcI+7r3X3fwELgQG5bZ1k0QBgobu/6+7rgHsI77kUIXd/DljZZHVBT/qkgB+tXsCSRvffT6wrZqeb2bzE19+C+nqbhlJ8f5ty4Ekzm52Yv6LUtHbSp7yU6YxXRau5iV/c/WtloItVCxPg/B64lBAELgV+A/wkd62TGOyfmMxoa2Cmmb2R6AmXnBYmfcpLCvgpNDfxSzOWAts2ut87sa5gtfZ1MLM/AA9nuTlxK7r3t60aTWa0zMweIKS5Singt2rSp3yllE60ZgDjzayTmW0H9AVeirlNWZP4g29wOOFkdjF7GehrZtuZ2TcIJ+hnxNymnDGzzcxs84bbwEEU/3veVIuTPuUz9fDTYGaHA78FyoFHzGyOux/s7gvM7F7CjF51wE/dfUOcbc2y/zazPQkpnfeAk+NtTna5e52ZnQ48AbQHbnP3BTE3K5e2AR4wMwix4y53fzzeJmWPmd0NHAh0N7P3gQsJkzzda2YnEMq3j4uvhW2n0goiIiVCKR0RkRKhgC8iUiIU8EVESoQCvohIiVDAFxEpEQr4IiIlQgFfRKRE/H9UizkGtDpewgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(x.data.numpy(), y.data.numpy(), \"ro\")\n",
    "plt.plot(x_.data.numpy(), net(f_).data.numpy(), \"-b\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 110,
   "metadata": {},
   "outputs": [],
   "source": [
    "from swag.posteriors.ffg_vi_model import VIFFGModel\n",
    "from swag.posteriors.vi_model import ELBO"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 111,
   "metadata": {},
   "outputs": [],
   "source": [
    "init_sigma = .01\n",
    "prior_var = 100\n",
    "noise_var = 0.0005\n",
    "criterion = losses.GaussianLikelihood(noise_var=noise_var)\n",
    "\n",
    "vi_model = VIFFGModel(\n",
    "    init_inv_softplus_sigma=math.log(math.exp(init_sigma) - 1.0),\n",
    "    prior_log_sigma=math.log(np.sqrt(prior_var)),\n",
    "    base=RegNet,\n",
    "    dimensions=[10, 10, 10], input_dim=2\n",
    ")\n",
    "\n",
    "elbo = ELBO(criterion, len(y))\n",
    "optimizer = torch.optim.Adam(vi_model.parameters(), lr=1e-6)\n",
    "\n",
    "vi_model.mu.data = torch.cat([p.data.reshape(-1) for p in net.parameters()])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 112,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0 tensor(263.2856, grad_fn=<AddBackward0>)\n",
      "100 tensor(246.7050, grad_fn=<AddBackward0>)\n",
      "200 tensor(233.9647, grad_fn=<AddBackward0>)\n",
      "300 tensor(283.3163, grad_fn=<AddBackward0>)\n",
      "400 tensor(230.0835, grad_fn=<AddBackward0>)\n",
      "500 tensor(230.9272, grad_fn=<AddBackward0>)\n",
      "600 tensor(227.8602, grad_fn=<AddBackward0>)\n",
      "700 tensor(233.2167, grad_fn=<AddBackward0>)\n",
      "800 tensor(234.0340, grad_fn=<AddBackward0>)\n",
      "900 tensor(235.4178, grad_fn=<AddBackward0>)\n",
      "1000 tensor(262.3733, grad_fn=<AddBackward0>)\n",
      "1100 tensor(235.5795, grad_fn=<AddBackward0>)\n",
      "1200 tensor(236.4392, grad_fn=<AddBackward0>)\n",
      "1300 tensor(273.7194, grad_fn=<AddBackward0>)\n",
      "1400 tensor(305.7097, grad_fn=<AddBackward0>)\n",
      "1500 tensor(240.2477, grad_fn=<AddBackward0>)\n",
      "1600 tensor(253.6443, grad_fn=<AddBackward0>)\n",
      "1700 tensor(263.3564, grad_fn=<AddBackward0>)\n",
      "1800 tensor(239.5156, grad_fn=<AddBackward0>)\n",
      "1900 tensor(234.6884, grad_fn=<AddBackward0>)\n",
      "2000 tensor(300.3623, grad_fn=<AddBackward0>)\n",
      "2100 tensor(248.0540, grad_fn=<AddBackward0>)\n",
      "2200 tensor(278.4770, grad_fn=<AddBackward0>)\n",
      "2300 tensor(231.6956, grad_fn=<AddBackward0>)\n",
      "2400 tensor(255.3054, grad_fn=<AddBackward0>)\n",
      "2500 tensor(235.3830, grad_fn=<AddBackward0>)\n",
      "2600 tensor(280.1711, grad_fn=<AddBackward0>)\n",
      "2700 tensor(250.0381, grad_fn=<AddBackward0>)\n",
      "2800 tensor(284.3923, grad_fn=<AddBackward0>)\n",
      "2900 tensor(242.1906, grad_fn=<AddBackward0>)\n",
      "3000 tensor(294.9906, grad_fn=<AddBackward0>)\n",
      "3100 tensor(241.4184, grad_fn=<AddBackward0>)\n",
      "3200 tensor(252.2780, grad_fn=<AddBackward0>)\n",
      "3300 tensor(250.4774, grad_fn=<AddBackward0>)\n",
      "3400 tensor(238.2356, grad_fn=<AddBackward0>)\n",
      "3500 tensor(282.3735, grad_fn=<AddBackward0>)\n",
      "3600 tensor(270.3519, grad_fn=<AddBackward0>)\n",
      "3700 tensor(246.3112, grad_fn=<AddBackward0>)\n",
      "3800 tensor(300.0522, grad_fn=<AddBackward0>)\n",
      "3900 tensor(231.6030, grad_fn=<AddBackward0>)\n",
      "4000 tensor(304.0490, grad_fn=<AddBackward0>)\n",
      "4100 tensor(237.9833, grad_fn=<AddBackward0>)\n",
      "4200 tensor(239.3564, grad_fn=<AddBackward0>)\n",
      "4300 tensor(274.0934, grad_fn=<AddBackward0>)\n",
      "4400 tensor(276.4989, grad_fn=<AddBackward0>)\n",
      "4500 tensor(250.7343, grad_fn=<AddBackward0>)\n",
      "4600 tensor(287.1213, grad_fn=<AddBackward0>)\n",
      "4700 tensor(233.9584, grad_fn=<AddBackward0>)\n",
      "4800 tensor(272.8073, grad_fn=<AddBackward0>)\n",
      "4900 tensor(250.9424, grad_fn=<AddBackward0>)\n",
      "5000 tensor(274.9344, grad_fn=<AddBackward0>)\n",
      "5100 tensor(233.2765, grad_fn=<AddBackward0>)\n",
      "5200 tensor(325.1034, grad_fn=<AddBackward0>)\n",
      "5300 tensor(275.2826, grad_fn=<AddBackward0>)\n",
      "5400 tensor(228.0100, grad_fn=<AddBackward0>)\n",
      "5500 tensor(242.4404, grad_fn=<AddBackward0>)\n",
      "5600 tensor(240.3985, grad_fn=<AddBackward0>)\n",
      "5700 tensor(236.6899, grad_fn=<AddBackward0>)\n",
      "5800 tensor(255.3155, grad_fn=<AddBackward0>)\n",
      "5900 tensor(237.9981, grad_fn=<AddBackward0>)\n",
      "6000 tensor(290.1969, grad_fn=<AddBackward0>)\n",
      "6100 tensor(234.5439, grad_fn=<AddBackward0>)\n",
      "6200 tensor(230.6493, grad_fn=<AddBackward0>)\n",
      "6300 tensor(296.7185, grad_fn=<AddBackward0>)\n",
      "6400 tensor(251.7043, grad_fn=<AddBackward0>)\n",
      "6500 tensor(235.2150, grad_fn=<AddBackward0>)\n",
      "6600 tensor(246.0012, grad_fn=<AddBackward0>)\n",
      "6700 tensor(245.2428, grad_fn=<AddBackward0>)\n",
      "6800 tensor(313.8044, grad_fn=<AddBackward0>)\n",
      "6900 tensor(256.0544, grad_fn=<AddBackward0>)\n",
      "7000 tensor(260.7643, grad_fn=<AddBackward0>)\n",
      "7100 tensor(383.2936, grad_fn=<AddBackward0>)\n",
      "7200 tensor(237.6385, grad_fn=<AddBackward0>)\n",
      "7300 tensor(226.9409, grad_fn=<AddBackward0>)\n",
      "7400 tensor(235.5924, grad_fn=<AddBackward0>)\n",
      "7500 tensor(257.4045, grad_fn=<AddBackward0>)\n",
      "7600 tensor(253.8358, grad_fn=<AddBackward0>)\n",
      "7700 tensor(272.3237, grad_fn=<AddBackward0>)\n",
      "7800 tensor(237.5028, grad_fn=<AddBackward0>)\n",
      "7900 tensor(294.0803, grad_fn=<AddBackward0>)\n",
      "8000 tensor(270.0817, grad_fn=<AddBackward0>)\n",
      "8100 tensor(289.0949, grad_fn=<AddBackward0>)\n",
      "8200 tensor(241.1746, grad_fn=<AddBackward0>)\n",
      "8300 tensor(227.1842, grad_fn=<AddBackward0>)\n",
      "8400 tensor(313.6406, grad_fn=<AddBackward0>)\n",
      "8500 tensor(329.8138, grad_fn=<AddBackward0>)\n",
      "8600 tensor(258.4602, grad_fn=<AddBackward0>)\n",
      "8700 tensor(234.9744, grad_fn=<AddBackward0>)\n",
      "8800 tensor(275.8043, grad_fn=<AddBackward0>)\n",
      "8900 tensor(257.9002, grad_fn=<AddBackward0>)\n",
      "9000 tensor(256.0947, grad_fn=<AddBackward0>)\n",
      "9100 tensor(264.6646, grad_fn=<AddBackward0>)\n",
      "9200 tensor(243.6844, grad_fn=<AddBackward0>)\n",
      "9300 tensor(227.9523, grad_fn=<AddBackward0>)\n",
      "9400 tensor(332.9117, grad_fn=<AddBackward0>)\n",
      "9500 tensor(275.6416, grad_fn=<AddBackward0>)\n",
      "9600 tensor(247.4072, grad_fn=<AddBackward0>)\n",
      "9700 tensor(263.2431, grad_fn=<AddBackward0>)\n",
      "9800 tensor(241.5271, grad_fn=<AddBackward0>)\n",
      "9900 tensor(239.6415, grad_fn=<AddBackward0>)\n"
     ]
    }
   ],
   "source": [
    "#     train_res = utils.train_epoch(loader, vi_model, elbo, optimizer, regression=True, cuda=False)\n",
    "#     sigma = torch.nn.functional.softplus(vi_model.inv_softplus_sigma.detach().cpu())\n",
    "\n",
    "n_steps = 10000\n",
    "\n",
    "for epoch in range(n_steps):\n",
    "    optimizer.zero_grad()\n",
    "    loss = elbo(vi_model, f, y)[0]\n",
    "    loss.backward()\n",
    "    optimizer.step()\n",
    "    if epoch % 100 == 0:\n",
    "        print(epoch, loss)\n",
    "    if epoch == 5000:\n",
    "        utils.adjust_learning_rate(optimizer, 1e-8)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 113,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_preds = np.vstack([vi_model(f_)[None, :].data.numpy() for _ in range(100)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 114,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-1, 4)"
      ]
     },
     "execution_count": 114,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAD8CAYAAABq6S8VAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO2dd3hb5dnG71db8t5xlhNCAgRCQwmbllH23qMpqxQohbKhlJVSCKu0UMomhQINM5QNH3tDIWETRgkj27HjbUuy1vv9cfvNOXYkT8m27Od3XbqsdY6ObOs+j+73GUprDUEQBCF7cQz1AQiCIAgDQ4RcEAQhyxEhFwRByHJEyAVBELIcEXJBEIQsR4RcEAQhy0mbkCulnEqpj5VSz6Rrn4IgCELPpDMiPxPAV2ncnyAIgtAL0iLkSqnxAPYFMC8d+xMEQRB6jytN+7kRwAUA8lI9QSl1MoCTASAnJ2fLjTfeOE0vLQiCMDr48MMP12qty7reP2AhV0rtB6BGa/2hUmrnVM/TWt8J4E4AmDVrll60aNFAX1oQBGFUoZRamuz+dFgrOwA4QCn1I4CHAOyqlPp3GvYrCIIg9IIBC7nW+o9a6/Fa60kAjgLwqtb6VwM+siS0twPBYCb2LAiCkL1kVR650wnE4xR0QRAEgaRrsRMAoLV+HcDr6dynHZcLcLuBSITXnc5MvZIgCEL2kFUROQB4vYDDAYRCgLRSFwRByEIhVwrw+Sji4fBQH40gCMLQk3VCDtBS8XqBWAyIRof6aARBEIaWrBRyAPB4KOjt7UAiMdRHIwiCMHRkrZADgN/Pn6HQ0B6HIAjCUJJVQh6LdU49NH55IiEpiYIgjF6ySsjjcUbfdl/c5aLNEolQ6AVBEEYbWSXkWjPyDoU6++JeL/3ycFj8ckEQRh9ZJeSJBKPuYHD9Un2fjz8lJVEQhNFG1gl5NEoxD4U6++IOB8VcSvgFQRhtZJWQR6MUabPo2dra2Re3l/BnvV8+fz4waRLPUJMm8bYgCEISskrIc3OBnBzaJ7EYf7a1dfbFfT5qXzicxSX88+cDJ58MLF3KN7F0KW+LmAuCkISsEnKtubDp81lRd1sbL3bRzvr88osvXn8RIBjk/YIgCF3IKiGPROiBuzp6Nhq/vLW1s2hnvV++bFnf7hcEYVSTVUIeCPDidlteeCRCQW9r6yzaWZ1fPnFi3+4XBGFUk1VCHo8z2vb7gbw8CnVbm7UI2trauVgoa/PL587lGctOIMD7BUEQupBVQt7SAqxda3nlhYW0UNraKNahEMU8Hre2MfnlWdW/fPZs4M47gaoq9iGoquLt2bOH+sgEQRiGZJWQ19QAzzwDNDUx8vZ4gNJSirrxyUOhzoufxi/Pun4ss2cDP/7IA//xRxFxQRBSktZRb5mksRHYbTdg9WrePuQQ+uQuF1BSAtTXM7HD+OFOJ90Ipfgcr5dC7nDwBCAIgjBSyJqIvLAQ+PWvGaCefTbwwguMuiMR/iwupnfe3k7Rb2rqXK7v8VDQ29s7Wy+CIAjZTtYIOQCcey5w0kmMun/7W+C11yjQJs2wsJCLoJEI0NDAi91OMcVCWeWXC4Ig9EBWCXk4TAH/5S8p1iecALz1Fis+jUAHAkBBAcW+sZGLo5EIt1dqBBQLCYIgdCGrhDwaZZR97rnAQQdRjH/5S+C994CiIkbn4TAtlPx8btPUxEVSk5ZoLxaSTomCIIwEskrIm5uZgtjaymr1vfZihsrhhwOLFnHR0+9nNO50Wr54MAisWWPZLKZYKBqV4c2CIGQ/WSXkVVWMvBsbKdZz5gA778yo+8ADgc8/56Kn18uI2+ejoJsmWnV1lph7vRT0cFgWPwVByG6ySsjr65ly6PNRlN1u4PLLgW235WP77Qd89x3F3OFgtG0KgkyFZ0ODFYXbFz+zqvJTEATBRlYJOUBrxelk1kl9PUX76quBmTOB6mpgzz2Za15ayudHo7RbjHiHw1ZE33XxUzJZBEHIRrJKyN1u9iPXmh53UxNFe+xY4C9/ATbZBFi+HNhjD/rpJSUUbK25nck5DwYp5qZ3i6n8lMVPQRCykawS8rY2LnQ6HEwzzMtjeuGSJcCUKcCNNwIbbMDbe+xB4S4pYbStFNMUIxGKdlsbxT6RsCo/zeQhQRCEbCKrhNznoyA3NdHrzs2ltVJfDyxeTDG/5RZgwgTgiy+AffahDZOfT+E2/rqp7mxp4cVE+KY1rmSyZDEyIk8YhWSVkCcSFOWCAkbmq1ZR2EtKKO5ffcWI/LbbgIoKYOFCLoDm5FDAm5ut62ZIRVMT9wVYWS6SyZKlyIg8YZSSVULu9TJqDgQ4YyGRoLUSjzM6b2gAvv0WmDyZYl5czMrPww7j4qfbTeHOz+8cfTc1WZPV/H7JZMkKkkXeMiJPGKVklZAHg4yqlWLkPGECLRHjdft8FPMff6SY33wzo/cXXmAX2PJybtvYyCje6bR88cZG/pRMliwgVeS9dGny58uIPGGEk1VCboYtm8VLjwcYM4aCbBZBnU7mmK9ZQ8/8hhtopzz+OHDKKRTzeNzKanE4uN9QiF57JGJNIUokpCfLsCRV5O10Jn++jMgTRjhZJeQ5OdbiZnMzBbigwBqk09hoReu1tbRMpk8Hrr2Wtsy//w2cdx7FPBzmZ7+khNuYTJb6eqvEX3qyDFNSRdjxuIzIE0YlAxZypZRPKfWBUupTpdRipdTl6TiwrrS1AXvvDbz9NjBuHG83NFC8c3OBzTazcsvjcQpzbS2ft+WWwBVX0Be/4w5Wg5aVWfM+S0v5fBPx19fzttvNE4CZCSoME1JF2GYknozIE0YZ6YjI2wHsqrX+CYCZAPZSSm2bhv124o9/BD79FLjwQuD8860pP+EwuxtqDcyaxV4sLS0U5dZWFgyFQsDPfw5cdhkj7RtuAK6/ns9taqJoFxdzH7EYo/26Ot4vaYnDkO6GU8uIPGEUMmAh16QjgQ/ujktalwhfew34xz+s22+9BRx5JPDuuwy83G5G5q2twOabcxE0GuVnuamJ1Z6trcAvfgFcdBE98KuuAm69lRksDQ3cR0kJtzE55g0NFHefz2qwZUbJCUOIDKcWhE4onYa0DKWUE8CHADYEcIvW+g9JnnMygJMBYOLEiVsuTZVh0IWWForzjz8mf3ybbYAzz2T+eCRiFf2sWMHUxPZ2im8gwFL+/Hzg6aeBa66hSF9xBXD88RTpsjL+rK/nvr1eRu2FhbwdClk2bKp1NUEQhEyhlPpQaz2r6/1pWezUWse11jMBjAewtVJqsyTPuVNrPUtrPausrKxP+99119SPvf8+s1EWLKDAulyMxisqaJfYZ3XW1/Oy117ABRcwmLv0UuDeeynaa9dyH0VFjMzNyDiziCo55oIgDEfSmrWitW4E8BqAvdK1z7w84J//BJ59lqmGyWhpoef9+9/TL3e7eRk7ltG000khbm2lCLe0sBfLuedy+0suAe65x8p2yc2lmNtzzJuaOueYB4OSYy4IwvAgHVkrZUqpwo7rfgC7A/h6oPvtyj77sOT+kENSP+e992iTPv88RTYQYNFfcTEfN/3Ig0EK9B57AOecw8fmzAHmzaN1snYt0xqLi3k7GLRSHk3DLkDEXBCE4UE6IvJKAK8ppT4DsBDAS1rrZ9Kw3/XIz+ei57x5XJhMRmMjFzQvu4yC7PEAU6cClZW0WFparH7kWvMEcd553HbuXK6ZmWlC+fmWzRIM8j5TeGQvGBIxH+VIoy5hiElH1spnWusttNaba60301r/OR0Hlgxjmey7L/DKKxwikYpnnuFg5oULGVVXVfESCFDIV6ywslDsnvn117NPS9uDT6J5461QMD4HRXttg8STTyEY5MmhrY02jN/fTcHQYH64RUiGDmnUJQwD0pK10ldmzZqlFy1a1OftTHMrkw4YCnGh8soraXskw+nk5+r44ynU9fVMR6yrszJZcnP52Ouvswo0kQBOcvwTZyeuQRGakItWNHrHouGyv0EdcAByclgdmpPDhdVwmMdj/PN1H257GXkgkJkUucF8LWF9Jk1K3uOlqip1qpUg9JNUWStZJeShEMU8FrOmBQEcunzOOcB//5t62y22AP78ZwqwEfOaGor4mDFcVAWYm37t3BhicGE27sPFuAIlaEAAQTSOnYGG/3sfWvP5FRXUzEiEnrtJfRzUD7cIydDicCT31kzfB0FIIxlNPxwsTKWl32+lBkYiwKabAg89xF5KXm/ybT/+GDj6aFoyxcUs8y8tZTRdXW1Vg26zDTAHF8ONCObjWFyEa1CLYoTgR8GqxSgs5Ge0pYWNuYJBHpMp5Q+HkboXSCa68A3ma41UBmJNpWoXII26hEEkq4TczNpsb6dwOp28LxRihHzSScAjj3B2ZzJaW62FUL+fQWtxMe+vrubPWAzYqHAtLseF8CGIJ3AozsLNWIMyRMZNQWEhUoq5x9PRl2XChskPIBMf7tEgJJlcAxiox91duwBBGCSySsidTkbl9fVWNK41hT0apU0yYwZw333ACSdQbJPxwguMzr/9luX85eU8QaxZw4jaccwvUYVVuBx/RB6a8DL2wG/UvVhx6pWIRq3URKXozRsx93op5pFLr0C7v7DziypFkUi3EI10Icn0YuJAh1FIuwBhGJBVHnlzMz1yh4PibUayORz0p3NzqWFtbSzsee89djpcsyb5/pQCjjoKOPhgnhyqq7mPigrA9e6raJ//GOqbE5ir5qJOF2OzzYC772Y6o8PBiLy+nsdQUGB55uEwEP33w/BecQk8y5bwhey/53QvRprpOMuWMRI3zaNGAulYA+ju9yMet5BFjIjFzpoaaxqQ08nPYCJB4TQ9VfLzGWEnEnz+6tXAddcBL76Yer8bbQScfTY/z8uXc8Fy0iRmojQ10bq5/nrOCK2qYsS/+eZWtajpYV5QwIXTdWIeBbybTaWYd0UWI3vHQIW2p6weWSwWsogRsdgZCDBTxWSrmAHKXq/VyKq6GvjyS4q4GTpx9dVMUTTbdeWbb4CzzuLPKVMowN9+S5HOy6Own3cesOGG/MwfeiizWwCeOEpKrA6Mq1dzO5+P97UvW4MI3Ou/qCxG9o6BrgH0ZJ2MdGtKGBVklZBHo7w4nexUWF5OITXReSDA+5Viwc+331pWzN57cyF05szk+w6HgRtvZNXouHHc14oVtE/8foryaacBP/kJi4KOOIL9X0wqYkkJTyjGomlp6RDziZVoh3d9MR9Ji5GZZKBC21NWz2B43FKwJWSYrBLy3FxG2VrT7nA4GGWXlVndDnNzebukhM9rarJmeDocwE03Ab/7HW2TZLz3Hqs8Gxq4v9paRv55eRTqE2a8gx1c/0VbG3Dc7Cj+eeoHUMoaQxcIcNvqar6276rL4Pa7O4u5RHy9Z6BC25uIPpPDKKTyUxgEssojN9kpACNorZkl4nBYszq1tsrmTddCU4rf3Mxv1W43o+0bbkg9eB0AfvYz9mJJLHwfY168F77WpWhCKRSA53AAnsahAIBz9/oCf3p0M2htjaBrabEWQAuemY/2i/+M6LLV8EwcA+9Vc0bOYuRg09eF3aGufBUPXkgjI2Kxs7mZYm6SQCIRqzTe7aaAa83HXS7+DAYprG1t3MfatYzQTQfExx5jOmIqinPCOKztnxiPJShHDYrQgCYUQEHhfeyAeTgNAHDYYcBdd1m57Q0NPInk5bFhV2GhdSIyBURCH+mvKPdF/NOdASRZMUIaGRFC/vHHFEOPh/6zz8fPiInM8/P5uTGRuRF1I/rt7fzZ1MRLfT014ZNPqAVNTalfe1u8gZ3wEiqwBhVYjVYUwAFgKTbAtbgEUXix/fbAA4f9B0V/uxThFWvQMGY6mk75A/yH7ovKSlov0fseROSyK+Fe/j18VRUjK1UwQ9TWcnH589lX49O2KlRjLHLQhgI0ohCNqPLWYYOcamxQvxAbjIui5JrzoX7Vj99pqhPFcccBzz3XP3GXiFxIIyNCyHfbDfjiCw5S3nFH9k8xYg7wM5efz2jc6bQWQU1QZKJ5e6QeCtF6Wb2a/vkHH6R+/RLUYl88ho3wFcZhBZqRBzcSaM2djCscl6O5GZiG/+FRHIQqrEAMTtS7x6D+ousROHRflL+5AGXnH49oKIoIPHAjCl/AKQUkSaivZ0O0228H/ve/vm3rRRiVpVFsslUeNt0U2GwzFopNn97RCycVqUS3P3UAJrJfujTzdQTCqCHrhTweZ2qg/XPm8wE77ADssguw7baMygF+TkyHRFNtaao8IxErD11rpgpGo/TQYzHgqaeYvRIKpT6WrfEm9sDTmIilaEMp9NG/gmu77XHNWauxKlGJUtTiQRyOrfAR4nCgfsymaHjqHTj32xOVNZ+gHGsRhxPt8MKFGPwlOVylHYkFPX0gHqdoX3ghA+B0D7pWChg/nr15pkxhYde0acDGG/PX7nSnsEGS0V1EnSyyN2JeVTVq/77CwMl6IQcowgsXUmyff55dDw1OJ7DddhT1nXbiQqMRb6WszoQ+H4U8HqdQxONcmIzHrYXT2lqKyccfpz6WYtTiAPezmHXEJES23BmRCJB74a9xK87El/gJfAjhdpyEA8AZGw3fNaJ+yhbQ0KhALcagBgpAGD44EYcfIazrKDBKIjbzLam6mn/Xa68FPvpoaAZ1OJ3AJHyHWfH/YissxNb4ED/FR8hBMPkG3XncYqcIGWJECLmpsgwEKLorVgBPPgk8/TTwzjtWBOdyAdtvD+y3H7D77ozKzfzNWMyK1H0+foBjMS6AtrXx82kee+AB4NZbuV0ylAK22go48ECeKGKXzUFJ+Hs8gOPwBnaDQgIX4CpcMO5hOL74HE3Tt8Pa1UFE4UYFalCJNXAigTB8cCCBAIKWmI/QD735Xa9eDfzwA7tWvvBC9+sTQ4VCHFPwHX6GN/E73IZZ+Mh6sLu/jyxwChliRAj5998z68SkGvp8lo0SDLJA57HHgLfesj4vgQC99UMPZYtap5OC3d7O57hcfI7XS8+8sZFWixnntnw5+5h359OWlnKW6OS6/yL02DMoxwq8gV/gURwDANjnJyvwj/+MR8nLD6P59xejPuJFGwKowFqMQTW8iCIMHxQ0/AjBAT2iPvShEP9uy5dT+557Dnj5ZVbf9vbfz+MBJk+mz73BBvz7RCJA80f/w8q3vsNyjEUNKtGIIsSTVdKmgQPxH8zFJdg0sLT7b0wSkQsZYkQIeThMwW5r4/VQiKKrtdU4y++n7/3888CCBfyqbqioYJR+5JH0RE3bWqeTQu528z6T5miKhhIJLrw9+mhq31YpLsL+ovhdRJ5/A4XhH7HStzlujZ2C9pgL06cD//oXsOFHjyD05+tQv6oNTWVTURReg/EtixFACGFwxJAfITirJmT1hz4SoUW1ahUj7w8+4Pi9JUtoY/WWadOYNLLXXqzkzc3l38rUDigF6PkPQF82B/HlK9A+dgpWnn4Vvp56ABYvZtuFb75hlW9Ly8Dfl0Ich2y1An+4pQpbbZXiSUOduy6MWEaEkBu0prgmEhSMcJgf0nCYt2MxftA9HkaBTz9NX335cvsxMErfcUfejscZ2fv9FIfaWms2p9vN/X3+OfC3v3XeT1cqKphTXlbW0dI2Atx8MyPS4mJmYeyyC1+vvh5oePA55Fx9KSbEvkM+WtAOH7Q/B77bb4Tr2F/2+3c0FCQSVr+ZFSv4DerVVznMo6Ghb/v66U+Z9LHppszDz8tL3Za4J0z6afUG2+GzVYX4AVOwBBtgCaZgCaZiOaoQRopGPN2w3XbAb38LHHQQs6U6MZI7UgpDxogS8lRobXnhra2WhRKJ8PLhh/xa/8orVlZKfj7b2O6/PwUDsIY8r11rLYQGAhSScBi45x7uJ1Vk6XBwwXX77Rk95uYyGPv6a0b5l1wCnHKK9Rp1DzwL9+03Y0LNhyiZkIfon+YifvhR6xpvDXfCYfreK1fy8vXXtE7eeSf1+kIynE7aXxdcwEh83LiBCfh6JPGuNYA4FFZhLD7GFlh0ydN4/30udK9d2/vj3mUXNl7bZ580Hq8gdGFECLmJwu3ZKPbrqbaJRvkt1xQC1dVxge3ZZ4HFi63nbrEF+5PvsYe1ENrURHFqbaUImxTHjz5ib/JVq1If75gxtHLGjmXvl6ee4kkEYBR33XW8v6GB3wAAPreigscdiw3vKtDWVvawWbOGv4dvvqGl9f77fUsdzMvj7+OEE6y+OSUl1N20ksq7NnTxsGtraafNmdN7Ud9yS6avmm96gpBORoSQmx7fqegq7MkupuFWOEwh+vxz4MEHKeqmjD8vD9h3X34T3mQTCsrq1fwwt7dbnnpNDbd95ZXU65IOB3Pdd9iBZfpffsltolH2Qb/jDhastLVRECMRClllJbc1WTY+3/CI9Ezu/apVPAGtWAF89x2/obz7bu8FXCkuXp5yCmsAysuBoiLaTxn7FpLMuzZ042GHQsAtt9AdaWzs3Uvttx/XVYqLB3jMgmBjRAi5KdoxkZoRNlOOb7/e9dKVaNTKXDGVnk88Adx/P/DZZ/ZjZUbKrrtaA59DISuFMRJhdH7PPd1H5+Xl/Nq9wQa0a+65hyeG/Hzgmmvo18fjFPOWForahAnWUGen0/LvhwIj4NXV/Eazdi2D22eeAV5/vWPodC9wOPg7PeMMWsfl5Yy+CwoGyUayV1w6nfyl97JIJxwG/v1v4KqruIDbE5WV/Ibyk5+k6diFUc+IEPJkEbnJXjCl+PbbdtFLJu7xOEXZnCDM/hYvBv75T6YyGo93/Hh66TvtxONoauL+TepiYyNz2l98sfvofOutgZ13ZoT96KOWtXPssfwKn5vLr/T19RTuiRMp9ub4/H7qz2CRSPDbQk0NT2Jr1/L4Xn2VVlF9fe/2oxTthlNP5e/SCHh+vmVXZQuJBPDGGxT1BQuY5ZQKt5vZSr/MrnVrYZgyIoTcZKsYITbXjXfe9a3YBT6Z2BtMQy27/+528wN6++38xm2i7cJCZqXsuisj1KYmvq7fz5PCypUcYLF6der3UVIC7Dn1S0z56FF8EtkUz+JgJODEZpsxw2WTTXhiqKvjMRnhM+9zMBZB43EKeGMjxbu+nr+Pzz5jodS33/Z+XzNmAKefzglL5eW8FBR0bp2QrcTj/CZ37rnd2++XXMJ6hGx/v8LQMiKE3GSgJPO+uy6MJRN7O12F3aSo2bsmGjsjHGb0/Pe/W2X7Ph/tkAMPtFIJXS5aNKEQffPXXksdnSsksDk+wM54DhHk4UGciEYUIzeXkfns2Xy/NTV8/fJyLoQ6ndxnphZBYzEKeFMTI++GBgr46tWMwF97rft1Cjvl5cCZZzKVsKyMi79FRUNrEWWKeBy47TZm3KTq0/OrXzE6H8xvVMLIYkQIeTxu9RzvGp2nqojuarPYvXSgcySvtdVAC+B29n7noRDw9tvMSjA9zB0OFqsceCBFKhLh/c3NzDe/7z5G6anIQTN2xP9hE3yJl5xHYHF8OgDggAM4a7SwkL55czMXYauqrH4x6fLNzUmstZX+fH09L42NfN1332XUWVfXu/05ncDhh/NSWspvFCUlPP60Z6IMM774ggvlqSbM7bwzF9a7Tq8ThN4wIoS8J1LZLvafXVlPWB56CLE5VyC0ci3iYydC/+FCOA47FE4nP3wmLfDrr4F//AP4z3+sCHWHHSjA48Zxv6YC9aWXuCjYXcV9FZZgFzyD2n3Pwgsv8DXGjwf+8hdOKjJpk04nF0GLiqxvDv31zeNxq1K2rY0ibmyUlhZmpDzyCDNtgPW7sSZjgw2Ac85hZ8HSUv4uCgp6aB87wmho4EzXl19O/vjMmWwjkZs7uMclZD+jQsh7QzJffd31hx5mOkUoCA2gHR5EvQVIXHcdcMhhUIqWhssFhOc/hvjV12LlKo07c87Fg5HDEIqypn/mTEboU6dS5E2/83vv7d5HdSGKLbdxY+ONGfFXV/OEcOKJwPnn8/qaNdynPUUR6L1vbqphjYCb67W1FHEzqu7552mjxON8jZ7avjgcXDv45S+5gDlhAoXcDPsYbUSjnA07b17yx6dPZ8FUYeHgHpeQ3YwaIe+ahmi/r+vj62274VTo5cuQgGPdJQoXgmOnIf7O+5Z3/uRj8Jx/FmLRKGJwIQYn6pyVuHfne3H/Bxut6+kxfTpTF6dPt6Lzt98G/rMgjmg8dQhdUMAovL6etgbAnPPrruNJoqaGEXNuLmtcvF5rvF2yfHPzLcK0MwgGrQrYNWso4KZ3zZtv0kZpbeW2JkOvO8rK6A1Pn8686aoqfmPI6ig8DSX2iQRw+eVc5EzGVMcSvH/rRyg65Yg0HLAwGhgRQm76hyfLEU/lk/eJ/HywaLszGgrtNc3r0h/1jjtCrV4BJ+IANCLwIg4HohUT0bTgFdx3H/Dww1Za2vTpTF2sqqJ4NjQAj9xRhy9WlnR7OJMmsdfIO+/Qr3a5WM9y1lk8jtpaRrumctSMuTNVqdEoL/G4lZkTi1G0jYBHInz8ww+ZjVJdzdc2WTg9sdVWwNlnU8ArK2mlFBZm+YJeb5pe9VLoEwlmPp1+evL/z2n4Cgvv/Az5Jx2ZoTcjjCRGhJB3l7WSqmzfHp12V86vFNYr4V73m5nI0m2Tdx7PL0ACCgoaDmi4EEEEXoThR/sP1UgkKNZ33cXPu4lup0+nhz5mDN/H118zV727XtwOB8Xc6eRsUYB9SObOBTbfnGLe2soofuxYa4HW4+nokR6zmomZdrJ1ddYQ6P/9jznzX3/NfefkWILfHQ4Hc98POYTbmOKetPZGGQgDiah7akPbx+6G8TizVU76TRwa65/hNnd+jveaZ8gCqNAjI0LIU+WL9+YtdBWXpAL/wHzglFOgQqzVV0DSD2ikairal1UjAQUNBScScCIO14SxiHy5BPG4dZwrVjA3fP5863M/YwYzXcrKKLJvvMECm+7eh89Hz33pUkb6Tif96AsvpChXV/N9VFRYVovWFHWTkbJ2LQXa6+Xzb78d+O9/uX+/n2+1N5kpBQV83RkzuN2UKXzdYWOlDLSNbE+DIfrRbzweB+51HYPf4F9JxXzrrfl/MGx+h8KwJGNCrpSaAOA+ABVgEHun1vrv3W2TqV4raeHhh2lqLl8OTO4ITL4AACAASURBVJgANecy4MgjO0f/Dz0I/bvT0B6OI8ExEIA3AOdNN8DzK/qdZpSc2WbpUuahP/igZVlsvDGHXkyYQAF9/PGeC23y8mijLF1KrTG52occYqUPut2MtkMhvrbxx+1FTq+8wmN0ufj6y5b1rk/4RhsBF11EK6WoiFkqZWVW7/ZhwUAHO/S0fT8nAMWrNsA9y3bASbgXwPorwHvuyVz9bKt0FQaPTAp5JYBKrfVHSqk8AB8COEhr/WWqbQZa2Wm9dk/HZm3Xn4t5raT3LVgAfdU1iK6sRqRyMvS55yKx/0HrCoyMR+x2W8VF8Tj14dZbmbZomnRNmsT84ilTKORPPNF92TdAQQesYQlVVezcuOmmfH2vl5ZHMEirZOVKpsOZgh6l+Hp1dd33CrdnrBxwAIc8uN30wydPpqAPCyvFzkBHrfUU0ff3RDF/PqK/ORU3h2fjHNwKYP1f3IknsrAoG9oXC4PPoFkrSqknAdystX4p1XOGY0Te1T/vzoqxY8S9vd0q8wco2qa4yOm0Bjs7HJaw/vOfjMKNaI8Zw2EF06axEddrr/XsVefk8LVM0yq3m02aysqoQ5EIfXC7ZTJ1Krf5/vvu9+1y8fUDAeCPf6TH73RSwDfccBgXtQw0Ige699gHYt3Mn4/2P/4Jly//Ja7G5UmfMmcOS/qH1bccYVgwKEKulJoE4E0Am2mtm7s8djKAkwFg4sSJWy7tLqE6BbFY/8dY9rQwmg4SCQqqyb12OGhvhEIUVNMx0URbbW30qv/1Lzo6pud1YSEFfepUdha0d2NMhalATTXOzOR2NzZ2X2kKWJONAAr2pZfyuAMB2kGTJg3ziHEwRq0NMD0xFGKe+b/+lfzx++8Hjj46y7N/hLSTcSFXSuUCeAPAXK31f7p77lAWBA0GpkWu1taYuGCQHnY4bKUIer1Wf/TqalZRPvSQNUrO7+fEnMmTmX/el0ZV5oThcvE1e9tmNieHJxil2EvmV7/iN4bycuawl5cPQyslGVkwaq2lhYVjr722/mNOJxfAd9xxdBZUCcnJqJArpdwAngHwgtb6bz09f6QLOWD1LzGRrakIDYUYFYfDWFcp6vdb5fJr17IXx8MPWymBDgczRMaOZdvbVH08BkJBAV8/GmWh0YUX0uJpb6eXPnMmj1NIH1qz6GubbTicoyt5ecCiRfxmlhUnTyHjZHKxUwG4F0C91vqs3mwzGoTcYPxzk8Hi9VpRcnMzo/Su3QzDYVZvvvceI/QPPrAspbIyBpirVnXfKre3eL3MgjFtemfMYDViLMaIftNNKeQSFWaGRIJrFVttlXz6UFUVsHAh/+6CkErI0/Hx3AHAMQB2VUp90nHZJw37zR7mz6dx7HDw5/z56x5yOKwcbYfDKpF3uWhTTJxIIdWawh4MWo2xDjqIjbnuuYfXCwtZAPThh7Rixo3j9v0R2bw8K9JbtYqveeKJ7OwYjzNCNz69iHjmMP8yTz6ZfN1h6VJ65d0VjQlCVhUEDUt6sbBmbx9gmmjF41Z/FK/X8tFbWmi/mCEXbjfWVYqaFMJ33wU+/bRzRkteHp/b1pZ8cr3TyROByZ4JBq20w402Yi76uHG0giorGZnn52fodyasRyjELKbf/z754+edx29KwzZTSBgURkRl53BDa0BPmozEsuXQUB3FQaz21BOqkPjy65TVmqb/ib33uVmcNDNEw2FrDJ3xSFtaGI1//z2n1X/+OdMLk4m32a/LZZ1M7OKfm8tIf7/9rP4sG25IK2VYZ6WMUJqbeUJNlsmiFHvbH3GEFAyNZkTI+4m9za39si7K7tJoi/1XElAAVCjYY18YM8zC+OimvN8MxDC9UkyfGdO1sLnZaj3b0sIF0NWrmS9eV8ev4o2N6wu818uFy1mzmBduUgkDAWbHjBkjVspQoTXXRvbdl/ZZVzweNlDbYgtJSxytpBJyKTnowAi2EdJUfV3WmwE6sQxq2Y/rxHsdVVVAL/pmmAVQk4poxDwWs4qf/H6rEVYsxtsFBbRColEKeW2t1Z7WlOa7XPzARyJWhaZJeywqsroUlpRYE3wkO2KADCDtUSn2cL/nHmCXXdbvexOJAHvvTTGXTBbBzqgUcq0twbaPjzMYkXa5Os/1TBqpXvXn5B753Ll9Pi7jiRtbw94WwOHoXC2aSFgtaisqGE0rxfuamxnFm5ORx0MbxeejGHg81n6Li3m4mZj/Oeroul6ydClvA70Wc6eTaxbz5tFG6VrJvHYtsMcezD0fN85qASGMbkaFtWLE2h5tG0xfFHuPlD5/MDJUfGI/0cRi1snGnGiczs7HbsdYNut6qGtrgdUMnzDXpRQ8TaSjNUAHLS3ADTewXD8Z06ax1mDyZGs6lFhiI59R45EbwbJfDHaxNgI4nKIZrYH4/Q8gdunliC9bicSEKuCyy6COOhJOp2WV9OcDa34vZspR16HUQhoYaLOuLtTXA6eeyorfZBQWMmbYYw9aY/ZqXvm7jkxGtEduItauwm2aVRnhHo7/3GYMWywGxB94GDjjbKhQG5yIw7P8f3CecSIc/tiAI3zjmQsZZOLE5BH5xIn92l1REXD99WzZ8N576z/e2AicdhqbpF11FTOOCgu5hmIsOlkUHR1k5ZcxM0A4FOJXUNPlD7BK3nNzLe93OEUodsvDTK83fVk8V1yKQGgtctEGP8JwIwZHqI22jTD8mTt3/UTvfq6XAPyfHTcOuOMOlvGn4tNPmUJ6wgls4bB6NWsEWlv5/2WsNWHkklVCHol0Fj+zkDfchdukDLa1WY2zYjF+E/f52KgqJwfwLl8CJ5J8Be9Pc5Vuqk2FDDF7NgvBqqr4D1hVNeCOiw4HsMkmHAZy0EGpn6c1C8V2243zQZcsYSZTU5Ml6GYBXBh5ZJVHbgYFD8QrHgw62SUdVo8ZLmG/rEe6FssGo42rMKjEYhTlv/8duPZa6xtoKlwuCv+ZZ7I2wAQ5pnmb8dKF7CKTvVYGDRN9D8cV+njciro72SUefohyc3nsxrNPSqqv5vvs07fo+uKLO4s4wNti0WQtLhcXNC++mBW9O+3U/fNjMWDBAkbo11/POKCxkf+bwSAvxnYRsp9hJofZham6bG21fHr7mLWcHF7v9YJTsq/mxx0H3HuvNaTT5CZ3J+aprJhM9L8dbQyxZeX1cnHzqaco1FOmdP/89nZ67HvvDdx9N1s7NDVxfSkc5s/WVv7vio+evWSVtTLUGMvEXkRkLBOXK0PefH/sljTmMws2hpllFY+zkdqCBcA119AT74nKSvaa32orFpLl5lrZLabfj8czfNaYhM6MmjzydGJyr43fbX5VJhfdiHdG6U9u8jATnBHDMD1BtrXx5efN45+4q6uWjE03BS64gCmLlZVsz2CvZBZBH56IkPeS7hYqh2SRdQAT24f7qLOsI80FP+kkkWDDre++A/76V1ov9pqKZChFD/3kk/lvVlpKQTc1F06nCPpwQ4S8G4x4R6PW59H0WjHiPWRIdD18GKYRuZ1gkINCFi8G/vxn4KOPet7G5wOOP5656OPGsf9Obq5VAWyyXETQh54RkbWSTswINnuWielE2K+FykzRdQG0pITpL8ccI/nhg02aC34yQSDABdDtt+fi5s03cxJVd4TDzFM/8UTOi12yhCeDYNAqYDOLopKLPjwZVUJuKkKNeEciVuMoU1BkJugMK2bPZsR3//38RNXV9T6DRUgfGSj4yQRKccbnRhsxyn7ySeC3v+15WMiaNcAll9CRW7iQNk1trfU5sVck90nQpTgt44x4a8VUVpqWr4Dld2ddC9As+GovDC8SCaYbGv/80kt7Z7e43Tw/HXQQF0PLyxnseDxWW+WuowqTItbgOkzr6YFYVKPOIze9u03Bg32U2rCLuHvLMF5sE4Y3kQgj7vp64IUXgCuuoFXSE+PHs+R/5kz650VFtB3Nv6L5LKUUdAk+EItZVemANSimP4wKITd+nlm0tC/UDLnXnQ7kQyEMAK2tEYE//shS/5df7t22++zDZZnKSl4KCy3hNtGl1kkEfZQGHyb6Ng3LTCA5UBdgRC92xuPWYoxZtDTNqMxQ4RFBbxbbxI8UUqAURwROnMh5rX/9K/C3v1GUe+K554Df/Q549VVaNMuWseS/a6ZXPN7FQ0/VwrefrX2HMyaQNO0PIhFqj99PLTJjFzMRO2etkGttLVwGg/wH8nj4CwsEstD/7g09LbYZP7Iv5fzCqMPjYVQ9cSKw//7AQw8xn7wnGhqY0njVVcA33zC7ZfVqIHT/o4htMgOJ3Hw4NtkIrkcfRDTaIehzrob2D+9Mn4Fi5uSarB6t+a3ErCnE41YbD6NV6SbrrJWuX1lM0UJPK/KjArFehD4SDjMJqraWRUTXX88e/z2Rl8fpRT+LvoTCv1+JysQPKEQz3IhD+QPw3Hoj9FFHIxoF1CMPw/3nS+FZtgSqamQUp9nn55pWHfZKb3s1uKlTaW9n8FlRwSi9P4wIj9z8IgCZgJKUUepHCgMjHqdNUlcHfPklF0J7k9kCANu73sdpsb+gCA0oQT3GYA0K0ALHhIlwfLUYHo8lakoB3scegHvORVlZcZxMvO3Tx8z7BPh4e7t1SSQsm6WkhJF6fxgRo97M6DapMEtBmkeNCaMDp5Pi4vPxs3XTTXTj5s3ruc3tu7Ft8Bnuxln4C7bFO2hBDspQj5Ll9fC3W5an3w9E7nsI4dPPRjTUAi8UnMb6A4almKfqtWRvjhePM7g0Q92jUSvgNI6B389vMB4P78uEe5BVEbnQA5KzKwyQSIQpitXVnDh05ZX0wXvDDngN5+J65KMR+WUFKHznOeTl8YuixwP4Np8Gx/IfEIEXGgoeROBFZFhZf4lEZ/EGrEHldrE2102div35Lpc1g8DrtSJ586U4E+mHWRWRCz1gxFqaZQn9xOOhh+v1UnAmTgSuvpqi3j0a72AXfI6ZONdxI7b/9baI1lK88vL4jNDyFriRjxy0QkEhAg9icMG/dPmQZV2YqNs+wN0ItPlpv8+kW5pt4nHedruZ/ePzUaTN/kxM5XRaYygzUcciEbkgCEkJBllEtGoV+7bcf39PVosGQM9hl124GFpayiKikhIgZ9sZaF9VgziccCEGP4JQUFATJsL33eK0Ww5G2rS2LibiNl63sUxSZZIYD9xE5Ea8zWPGjnI6O88pAKzFT2PFmNcxnnp/kIhcEIQ+EQiwstPnY6+WjTZiEVF9faotqE5KAa+9Bnz4IXPPd9yR6XclZ1yJ0jmnw9UeRCvy0Ix8OD0uOM66Ao5V/AYQCFhdF9fttYvo2QuQuv5Mdp89ejYibG8zYIp1zHUzbwCwqjLN9C9jmxhxttsw9sVPgI+ZQev2k0ZJCaPzdCJCLghCStxuy2rxejnI+eqrga+/Tr2NsSCamzm5aLvtgBNOACb99EBELlUYd9ulGLvyK7SPn4q2C/+M2AGHIJGwimjM1CLAimS7Mw7MY8aDttsidkx0bI7PPG4E3LyO6fYIWAkWpp21edzYMGY7kxhmt2Hs+3a7LWslE8NoxFoRBKFXtLUBy5ez1/lddwEvvthzlaKxHPx+4IgjgF13ZTHShAns0GgWA03EajJAjNcMWJGzPVK3Z9TarRO7kJrbdqE20baxN8yipX10o3lNsw8j2vYJYeZijsncts8wyMQgGrFWBEEYEDk5wAYbUHzPOIMLoffdx3S7VBgRD4U4Q/ytt4Cjjwam1b6O4kfnoaThG/jKi5H4/enQe+3fSbzdbgq6PSI3EXpX+yXV42Z/5rmmlbX9JGDEv+v+7MJsj6bNicB+QrD3nBkKRMgFQeg1Hg+zBb1eimxFBXDbbSwmSkUoRDFXCvj+e+CquRq7I4yjEYEXPjhrauC5/E+Ia4XWXfZbFyErxdezD3ixR99A52jY7m8bG8SkBJptjWDbtzMzCEyDPTMNqatID+faFRFyQRD6hMNBe8QMYikspJh/9x1gz1yxEwpxu8mTgaU/JPAi9sL72AYnYB5m4z6Ux9ag4M7z4Lpgv3WNp0zrV5MdYsQZ6Gyx2H/ahd7unZviHHsfdXtUncmI2uSeRyL8mZeX/qKgtLg3Sqm7lVI1Sqkv0rE/QRCGN0pRwGfNou997rnA1pNXI5mIGxIJ4IcfgEn4FhvgGzShCDfifByCp/EIDkfDiia03/cwnNM3QkGFB0XbbYz85x9e14rD66UImmlefj/vM9GzicZNZG1y4fPyeKwFBbwUFnIffn/nyDsd2KeQNTRwoMeaNfzG0trKbwiZ6JaRlsVOpdTPAbQCuE9rvVlPz5fFTkEYOUSjFOj3tvodHm7eBS/gECTghEIcGsmbIeWjDjviNSzCDqhBJQBgilqC09WtODjxMNxIQEFDA4gXlCB22ZVwHnYw/H5G50DnBVB7KqF9ITJT2Mv3TRaL+WkwDf08HsuyGejCZ8abZimlJgF4RoRcEEYf8TiwzDUBH2AmnsfueBAnIwIfXIgghuQdohyIYzc8jTLU4BkcjiYUAQAm4zuchDtwNB5CHtiKMeotQPS6v8F1+CHIy+PCa6Y9a+Op28v27amH9sjaLtqZHGYz5EKulDoZwMkAMHHixC2XJmvuJPSe+fOlFF8YViSqJmPVsiDexyy8iy1xJ85BKwrhQRgR+FJuN8nxA35zUC2++88iPIkjUI9SAEAOWnE4HsHheBhbYRHax01B+xsfIB63bJKBDks3Qm1azZqf5rq5AFbEb7JXzOjIVLMP7CcC+0kgP38Yj3qTiHwQkeZYwnBk/nzok05GTciH97A1PsNGuAUXoAZj4UY7EnAgjuSrfD4fcGj8IWwTfQGfYBbewi/wLTZe93gVfsSBeBy7PnE2Jk/m4qnXa4l519RAu8UCdO4NbrJZ7I2vTDaL3aKx54N3lxNu34ddtE0lqb1PC2CNyesPIuQjCRkgIQxX5s8HzjwTNXUJvI9Z+BbjcTMuxA+YBpcjgZw8B5qaUm8+E4twDO5CAK1Yhip8iO2xCNuui9IBlrhvtx0wYwawxRbAtGlWi1hTuGO/2H10I8ZdI2oj2OY59nRDI9JA534rxm7pKuRmG/uJwZT0K8XeMyLkggyQEIY/8+ej9oK/YOGqMizzbYLbyy7Fp8vLoFTqtvmGHLTgGNyJ7fBfAAm4EcdKTMHbPz0D762oQk1N5+eXlgKbbAJsvjmw2WYU9okTrSIeu1edLE0RWF+E7eJtv921XYBdsLteUuWeD6TXSkaFXCn1IICdAZQCWANgjtb6n6meL0I+QCQiF7KEujpg0SJg5Up2T3z9dd6/8cZc3rG7g12Z4fwcx8ZvRnkB4D/2cFQesRvcbuDbb4EPPgA++wz45BMkjfCdTua6jxkDlJfzZ3ExbY3CQvrUeXm8mFRG+9Qee9aLvcio63UzGai93Zrb2drKcXnNzfxppi81NrLh2MEHA7/6Vf9+nyNi1JvQgXjkQhZRX0/BXbECeOYZYMECRrKbbMKcaxYSJcfnA/bdl7nqxcXs0ZKb29l/Xr2a/V8WL2bl6IoVbL3bV2lzOBgp28dI2n1283r2NrimaKkvnHEG8Pe/920bgwj5SEOyVoQsor6eQrt6NVvc3n03RXBSSRM2q3sZL2NvhBFIuf24ccCBBwIbbsjWuhUV1qKlsTCMsJpioLo6Tjpa8+InqH3lEzS0utHkG4PmCZuhxV+BlhZGz6GQ1W62P5jKU6/Xiu7NhCBTwGQKkQoKeFLaYYf+vZYIuSAIQ0pdHcW8oYHDnW+8kfZDMWpxGO7Fm9gTX2NGyu2VArbZBthzT8Yu5eUUT7vPHY9zny4Xvei8N56C48or4EiEoKBZd+ryAZf/CYm99++UWmgyTeyj3OwLpPYKUzNQwrT37dqTxSxX2dvwmktu7jBOP+wLIuSCMDpZuxb48kuWsK9cCcw5aRlWYSK8COII3IsgcvACDkQrClLuw+MB9tiDbXHLy62+JQ6HVazT1EQhDZx2LPLrv4cbEXgQhRftcCAO15hKON54vZNIm4XLrgMq7CmH9u6KdpLJqD1rxv54eTmj9v4gQi4IwrBg7Vrgq69og4QOOQIXN1+AzzALDsSxFx7HNHyFdxx7YmFi6273k58PHHkkcOyxtC1MdaXpd1JfD4R32hlOROFDGA4Ajo6mXg4k4Pzw4/U6IRrRNvcZkk0hsmfA2AXeLqlde6g7HMDYsaxM7Q8i5IIgDBuqq4H//Q/AC8/BedUVmIs/4HkcBADYDm9gp31z0DJpFh5/nAuX3VFYCBxzDHDUUdbkeo+Hl6Ytd0H7mrXwoh0eRKDhgIaCo7ISjrff7lSlaU8v7GqVmBRGe2ph18ES9kg9WdqhuS8TMztFyAVBGBJWrQKWLAGcLz0L/y1/w7yGA3AHTkcCTkyfDuy/P6cILVoEPPEEFyS7Iz+faX1HHMFFRqcT8P7fk4hdfBna4wk4oOFHGMoTQHjOlYjvdcC6BU7T+9wsWtoja/voNnuBUTJfPNnF7N9c8vM7pzr2BRFyQRCGFVpzdNyPPzIq9vmA559nAlYoxAXNI45gyqFSwNPnvopXoz9HvIcxCl4vc7WPO46lFYknn0Tkb7cgUVsDX0Uxci88DYHZh64rpQ+Hrf7nppTfXvZvH2xhJ9lcUPsJIJmIK8WMm0DqBJ1uESEXBGHYoTVb4K5YYXUP/Owz4IIL6KWXlNAymTYNKDrzENShFA/hWHyA7aF7MU5hm204Wm7XXbkIGgrxdYqKuG+/32ovax+8HAyyyMdkwtibZRnbxm6P2DNVklV/2q/n5kpELgjCCCMeZ1FQdbUV/S5fDlx4IX10vx847DBgy/+cj6K2JXBBYxXGYAGOwQfYrleCXlbGPPQ99qCQAlbet7FTTP63SS00WTDmYnLNu6Ymmn7jph1Ad+PhtB5YX3IRckEQhi2RCP3yujoKaSLBYp25c1lA5HAA+2/+I7b/5O8oQi1y0AaNBFZgCl6cfg7e+GZ8ryosnU5gyy2B3XZjbxazOOr1Mhq3D1020bqZRGT3zu1Ns7qmKto7J9oXRs3tnBzJIxcEYYQSDlPMm5oodqEQxfH++4Hbb+dztt1gDQ5cdROKwsuRWxBA4vAjEdl2F4TDwHPPUfRDod69XmEhsP32wO67W6Lu91vWSCTCi7FLjCCbCNz0XLFH6QZ7Xzu7oJumYeZbQV8RIRcEYdjT0sJ+KaEQbY5gkOL33/8Cf/gDfespU5huWFnJ7BSleL/Px0KjV14BXnoJqK3t/euOHQvstBMvM2Zwv6aC094cy/RYsfvgXdMJja9uny5kF/kpU+jR9wcRckEQsoL6ejb3jEQ6i3lDA3DiiawILSwETjqJi6AFHUWg8Tgjea25zaJFwKuvsmFXX2SushLYcUfg5z9nv3PTRCsQsHx1UwlqnyJkLyxK1i7X2DGlpcO0jW1fESEXBKE7qqvZYCsep5URDFIg/X7g1FOBd97h7aOO4gJmIGA10TJCG43yZFBbC7z3HvB//0cPvi8UFHCIxc9/Dmy1lTX42d5zxQyENoMpDF0LhwD+9Pv7P89ThFwQhPSR4e6biQQj77o6K6I1nQ5LS4EbbgD+8Q8+d7vtgFNOYZQeDlPAPR5LLM0iaEkJ8OmnwJNPAu++nUAs0bfUEZcL+OlPgZ/9DNh2W0bukQgfs+eeG3H3+zsPtjCeem6u1R+mr4iQC4KQHgapH34sxvNEc3NnG8PlYm/y115jb++WFnrcF10ETJ/OxdKWFivNz+GguEej9KYrP34ajZffhFcTP8NL2AsfYlavUhi7Mm4cLZitt6YF4/PRqzfj5kymitvNx0ye/Pjx0mtFEIShZhAnVEUiLBZqabEySUxhjs8HrFkDnH46i4hcLuB3v+M5pr6e0bxZaDSLlaEQ4D/nRJQ3LoYDTjgRQRh+vIVd8CL2wkJsjQT67nsoBWy6KSP1bbdlFozDQWGPRHjsTiePZfPNeSLqDyLkgiCkh0GeGRsK0WYJh62UwPzXnwD+fhNU9SrEK8bhpk1ux12vTwVA6+PGGyns9fUUc2OBtLUBLYfPhgdtKEEdABcc0PAhhFyEEYIH7zl3xcvTz8Ab347rsb9LKlwu4Cc/YWXp1ltT5J1OCvtGG0nWiiAIQ80QzIxtbWX0HQoB4ceeQfTKa1AcWwknNNrhB9w+fHDcLTj70e3R1MRKzltvZfTb0GCdX9xuoHXPg1BTHwegkIdmuBFDHC44kIAHMSgk4K0ohePpZ/D550xnfPVVLsD2F7eb9ss22zBvfeedmeLYV0TIBUFID4M9M7ZjYbVl6VrUjJmJUFsMoZYINDTy0QI/2tCCAiQqxqP1Py/j7LM5nFkp4Pe/p/XS0sJo2OMBcl95HLGL/oS6uB9tyIUH7XCjHQk44UQCDsTRDh/iL769roWt280TyaJFwBtvcP99ndVpx+nkiebkk/u2XSoh72fFvyAIo5bZsynaVVVUy6qqzIr4yScDS5ciD22oqP4YOS2r4EY7FDRaEUAQefAhBP+aH5Cby0rQM8/kod10E5tmRaM810SjQNtuByPn+jmYOAaoRDU8iEBBwY8QnIgDcKCgNAdFRUw/dDr5jQCgTXLNNcDbbwPz5gG/+Q3niPaVeJwWS7qQiFwQhOFLEhunFQHUoQh1KEICDrgQgwcxxCsmwvfWS4jFGHm//z5w9tlATQ1tjL/+lfng9fX0sMvKAO8zCxA+52LUR32Iwg0X2gGnH7h0DtwH7buuPN8UADU1MYsmGrVyyfPz6cF//jnz1V9/nZ5+d3i9QGOjlZfeW8RaEQQh+0ixsNoGP+pRiGqUwwEg4IwDc69E2y8ORCBgdS9cswb44x9phwAcDXfRRVw4TSS46Fj48gJ4r7oM9auCaC6fBvz+dOh9D4BSFHylKNTxuJX/bbJgIhHaNuEwDzMQ4EkjGAQ+/piR+3vvre+v77QTBb+viJALgpB9pFpYLSlBOFCM+uUtWFUxEzjtdBQcxQi6oYGC5j4+UQAADs5JREFUahpThUK0W268kT55VRWvT57M24EAL8XFFOPWVkuko1Er1RGwvPZwmGIeCFj2S1sbo+yWFqybPBQIMGd87VoK+8KF/Hn88cDll/f91yFCLghC9tHDwmp7O4W7upoRdmkphbtu/vPIueN65Nd8i8SYCQifeT6+2vAgXHIJ8M03jNZ//Wvgt7+lzWHK6EtKGM2b3O94nMKuFMXc57OadDU00KaJxRipFxZS1F0ubt/WRlEPBq2FUacTmDmTUXt/qjtFyAVByE56aAcQjVJUV63qKOF/YwHC512KNdF8eBFEIZqRcAUQvPw6RPfYH7ffDtxzD4V/0iRGxttsQwEPh+l5e70UX1Ni39pK8bZPCPL7uY+WFvZzaWnhba+XUXggwOdrbeXAh8Nc5OxP6iEgQi4IwggmEuFC5IoVgHPv3VG25lOE4EMtyuBFGAGEgDFj0fYcUwqXLOEEou++4/YHH8zbY8ZQbHNyGDGbhVOvl5F+JGKV3huhNxOFwmEKflOTFdGb3uVut1UvVVmZ/pmdkn4oCELW4/HQ1pgwAXCtWYZGFMKPMIrQiBjcCMMDVK+B202RnToVeOIJujYuF/D448A++wAPP0wBbmuzomtjk5hRcACfYzoehkL0xuNxHsP48TwhFBZai66mhW3XDonpQoRcEIQRgcdDW6R0XA6ciKMdXhSjDvloQQJuJMrHQj/1NHz7/wKRmVsAv9gF5274BB58kL51XR2HPh9/PPD112yj29xsjXgLBhmJm4ZXJvMlJ4dWicNhzfX0+ynkeXlWV0T7NKF0I0IuCMKIweMBCq48H/neOKLwIIwAxmElijxhJHbeGYE/nQ9HzUoACTTVtCJ+8WXYavkC3Hsv0xLz8jiN6LDDgNtuo1WzYgVF3AxhVori7fXSN29uZlTu81mThUwf8pwcZs+Y6N0MaU43IuSCIIwovMcfjcLbrkbO2CIEEUDT2BkovvkKFL/zLHzRRpSgHsVohhMJrIwWovHaO1BQwPFxzz4L7LcfF1Dvuot556+8QrGuq2OU3tREoTbFQCbLpaGBHrkZ3OzxWM/Ly7ME3T4WLl3IYqcgCCOStjaK67pc8bF5CMMPNyLwIoIg/FiBSqxFKXLeeR1FRbQ9nE5OILrmGi6KAsAOOwBz5nCh0uSWFxdbY+YSCSvv3Cxy5uTwefG4NefT4aCYp3tCkETkgiCMSIx37fXSuw6N2whehBCDGzG4UIhmbIpvMHlMFEox6jaZJ9tsw6zHs8/mPt55B9h3X+D2260MlJoa2i5NTVYuuakqjcWs4iCzMGqKijIRkYuQC4IwYikosPK528+7CCFvKTwIIwIPWpED5Q9gyl9+h8mT+VxjlTQ2Mso+/HDgwQeBAw5gtH3LLcCeewJvvskTRTTKqs3qam5jKj6NJx6N8gTR0sLrpvgo3aRll0qpvZRS3yilliilLkzHPgVBENJBURE9a8ehhyB67V8RHjsNfoQQHbcBmv46D9EjZqOsjFZJIMAUxgkTKLiRCEX5D39gEdG0aSw8Ou004NhjGZHn5vK5xkc3o+ZMuqHpy9LWxosp308nA/bIlVJOAP8DsDuAFQAWAjhaa/1lqm3EIxcEYTBJJBg5m4k/Ph8j6rY2Wh+mtL6ujsJbWEjBXbmSFkpjI6PpggJgwQJ2CGhspIAfdRRwzjks74/H6ckDVg652231NDf77m/mSiY98q0BLNFaf6+1jgB4CMCBadivIAhCWnA4KLReL2+bIp+CAj5WX88UQrPg2dREsd9wQw50njiR2zU0APvvz8Kho4+mMD/wALsZ3norTwr5+Yzs43HaKaEQ0xebm62GXGl/f2nYxzgAy223V3Tc1wml1MlKqUVKqUW1tbVpeFlBEITe43TSPnG7KdahEMW1sJDC29LCS14exba+nsJcVMRK0KlTaaMoRavm1FOBRx7h2La2NuAvf6GgP/64Fb2bDBXTgKu/M0B7YtAWO7XWd2qtZ2mtZ5WVlQ3WywqCIKzD7WaHRJ+PkXhrKy0Sv5/CG4sxevZ4aJE0NVlZJ2PGMEIvKrLK78vKKODz5gEbb0z//IwzgL32Yi9yn4/7NamIHk9m3lc6hHwlgAm22+M77hMEQRh2uN0U40CA14NBijlg9RaPxxl5m4VLrRllFxYyqs/NZU8V44vPmMF0xauuorh/9hkzXo47DvjiC+6zqIjb9TeHvDvSIeQLAUxVSk1WSnkAHAXgqTTsVxAEISO43RTV3NyO1MR2K0XQ5HybxcraWkvoPR4Kud9Pn72sDBg3joLucgG7705r5bTTuI8XXmC64kUXAd9/z4h/WKYfaq1jAE4H8AKArwA8orVePND9CoIgZBKvt3N/lHjcSg90OKzSeqeTGS/GM/d4GF37fIzm/X6eECorab/k5gInnEBBP+ggCv7dd1PQ77jDGuScTqREXxCEUU0oZI12a2+niPv91tSgeJzZKrEYrRVTNBQOW/f7/dyXyYppbLQe++Yb4PrrgY8+4mNHHMGsl/4gJfqCIIxM5s/nqB+Hgz/nz+/T5iazxONhNK21Je6mP4pJUzSZLWZBtKiI25pBEub+sjKmLJaU0D+/+26K+bhxjNbTTQYaKgqCIAwSXWd6Ll3K20CncXDdYdIJg0GmJebn02IJh60JPy4X7zcibhZEPR5G6Y2NFHKlaMEUF1u9yHNzaafssw/wi18w8yXdSEQuCEL2cvHFnQczA7x98cV92o1SzGJxOCjmxjdvb2d0bp/FGYlYPreJxHNyGMG7XLzU13Nbn48ngKIiintpqWW/pBOJyAVByF6WLevb/d1gxNzYKmZ4RHMzI3Svl5G7WRQ15fYATwBm7JtJMWxttaL23Fxub7ZLNxKRC4KQvZja+d7e3wNGzF0ua9ByaalVIGRmd/r9vG4KhswwCcDqweLzUcgjEX5JGGifle4QIRcEIXuZO3f9kfSBAO8fAH4/I2iTzWI8b7MQ6vPxOS0tjNhNymJ+vuWTR6NWsyyn07Jp4vEBHVpSRMgFQcheZs9mK8KqKipoVRVv93Khszs8Hoq1yUbx+63CH9Pe1uvlY8GgJdBFRXxOYyOj9kjE2hcgo94EQRAGHdNgyyx4Op0U6GDQ6tcSi1lNshwOirqZHGQyXkwK40CQPHJBEIR+4HDYpgy10zLJy7P6suTk8DlGuE0aoln0TCQYna9dy20zcoyZ2a0gCMIQMMDioFSYXHOv1+ox7vNRrE2Zv8NhVXM6HBRxv5/bKmUNgx6u/cgFQRCGHlMctHQpw2JTHJQmMQcsrzuRoDC7XIzIc3IYoZuhFIkEhdw03zJiHgzy0NKNCLkgCCODNBUH9YTLRavFLsyBAC9FRZ2HUhjrxczuBCRrRRAEITVpLA7qCeObO50s5Y9GeTs3lwubsVjnHPPcXKYher2ZGS4hQi4IwsggzcVBPWGKh9xuphiGw7RSiosZhYfDzGgxNouJ4sUjFwRBSEWGioN6wuSTmzFxTif7kufkWE22AGvxUyo7BUEQUpHB4qCe6Fo8pDUwdixtluZmRucmkyUTi53SNEsQhJHD7NmDItzJMIugoZBVCVpRYfUxN/1XpGmWIAjCMMa+CBoKcRHUzPg0tyX9UBAEYZhjXwRtb6d3XlxsTRKSiFwQBCFLMIug0SjFPD+fOeSZKNMXj1wQBCFDmKpOs9hZWJiZCUESkQuCIGQQt5tWi9aMzKWyUxAEIQsxBUEOR2Y8crFWBEEQBgGT0ZKRfWdmt4IgCMJgIUIuCIKQ5YiQC4IgZDki5IIgCFmOCLkgCEKWI0IuCIKQ5YiQC4IgZDki5IIgCFmOCLkgCEKWI0IuCIKQ5YiQC4IgZDkDEnKl1OFKqcVKqYRSala6DkoQBEHoPQONyL8AcAiAN9NwLIIgCEI/GFD3Q631VwCgMtGXURAEQegVg9bGVil1MoCTO262KqW+GazX7ielANYO9UEMEfLeRy+j+f1nw3uvSnZnj0KulHoZwJgkD12stX6yt6+utb4TwJ29ff5Qo5RapLUelb6/vPfR+d6B0f3+s/m99yjkWuvdBuNABEEQhP4h6YeCIAhZzkDTDw9WSq0AsB2AZ5VSL6TnsIYFWWMDZQB576OX0fz+s/a9K631UB+DIAiCMADEWhEEQchyRMgFQRCyHBFyG921HFBK/VEptUQp9Y1Sas+hOsbBQin1J6XUSqXUJx2XfYb6mDKNUmqvjr/vEqXUhUN9PIOJUupHpdTnHX/rRUN9PJlGKXW3UqpGKfWF7b5ipdRLSqlvO34WDeUx9gUR8s4kbTmglJoO4CgAmwLYC8CtSinn4B/eoHOD1npmx+W5oT6YTNLx97wFwN4ApgM4uuPvPprYpeNvnZW51H3kX+Bn2c6FAF7RWk8F8ErH7axAhNyG1vorrXWyitMDATyktW7XWv8AYAmArQf36IQMszWAJVrr77XWEQAPgX93YQSitX4TQH2Xuw8EcG/H9XsBHDSoBzUARMh7xzgAy223V3TcN9I5XSn1WcfX0Kz5mtlPRuvf2KABvKiU+rCjncZopEJrvbrjejWAiqE8mL4waL1WhgvpajkwEujudwHgNgBXgB/wKwD8FcCvB+/ohEFmR631SqVUOYCXlFJfd0StoxKttVZKZU1u9qgT8n62HFgJYILt9viO+7Ka3v4ulFJ3AXgmw4cz1IzIv3Fv0Vqv7PhZo5R6HLSaRpuQr1FKVWqtVyulKgHUDPUB9RaxVnrHUwCOUkp5lVKTAUwF8MEQH1NG6fhHNhwMLgSPZBYCmKqUmqyU8oCL208N8TENCkqpHKVUnrkOYA+M/L93Mp4CcFzH9eMAZM039FEXkXeHUupgAP8AUAa2HPhEa72n1nqxUuoRAF8CiAE4TWsdH8pjHQSuU0rNBK2VHwGcMrSHk1m01jGl1OkAXgDgBHC31nrxEB/WYFEB4PGOuQIuAA9orf9vaA8psyilHgSwM4DSjjYjcwBcA+ARpdSJAJYCOGLojrBvSIm+IAhCliPWiiAIQpYjQi4IgpDliJALgiBkOSLkgiAIWY4IuSAIQpYjQi4IgpDliJALgiBkOf8P/+8BOzkfTywAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# all_preds = np.vstack([np.load(\"preds_de_{}.npy\".format(i)) for i in range(n_ens)])\n",
    "# all_preds = np.load(\"deep_ensembles_preds.npy\")[:]\n",
    "pred_mean = all_preds.mean(axis=0)\n",
    "pred_std = all_preds.std(axis=0)\n",
    "pred_upper = pred_mean + pred_std\n",
    "pred_lower = pred_mean - pred_std\n",
    "\n",
    "plt.plot(x.data.numpy(), y.data.numpy(), \"ro\")\n",
    "# plt.plot(x_.data.numpy(), y_.data.numpy(), \"--k\")\n",
    "plt.plot(x_.data.numpy(), all_preds[:, :, 0].T, \"-b\", alpha=0.05);\n",
    "plt.plot(x_.data.numpy(), pred_mean, \"-b\", lw=5)\n",
    "plt.plot(x_.data.numpy(), pred_lower, \"-b\", lw=2)\n",
    "plt.plot(x_.data.numpy(), pred_upper, \"-b\", lw=2)\n",
    "plt.ylim(-1, 4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 115,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.save(\"vi_preds\", all_preds)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "py37",
   "language": "python",
   "name": "py37"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
